3 Parametric survival models

1 Learning objectives

Note

By the end of this lecture, you should be able to:

  • Define a parametric survival model and explain the fundamental trade-off between parametric and non-parametric/semi-parametric models.
  • Understand and evaluate the key advantages of parametric models
  • Differentiate between common survival distributions and identify when to use them
  • Compare and contrast the two primary modelling frameworks: Proportional Hazards (PH) and Accelerated Failure Time (AFT)
  • Interpret model parameters from R output (flexsurv and survival packages), including converting raw coefficients into meaningful measures like Hazard Ratios (HR), Time Ratios (TR), and Survival Odds Ratios (SOR).
  • Validate model assumptions using graphical techniques
  • Execute formal model selection procedures, applying Likelihood Ratio Test (LRT) and AIC/BIC, as appropriate, to identify the most parsimonious and well-fitting distribution.
  • Implement advanced modelling techniques, such as fitting Generalised Gamma models to evaluate simpler nested distributions or allowing ancillary parameters to vary by treatment group.

2 What is a parametric survival model?

A parametric survival model is one in which the survival time, \(T\), is assumed to follow a specific, known probability distribution

\[T \sim f(t | \theta),\]

where \(f\) belongs to a known family (e.g., Exponential, Weibull).

While the family of the distribution is assumed, the exact distribution remains unknown until the vector of parameters \(\theta\) is defined.

We use observed data (including censored observations) to estimate these parameters — typically via maximum likelihood estimation (MLE) — to best fit the survival experience of the cohort.

Note

Unlike the Cox Proportional Hazards model (which is semi-parametric), parametric models require us to specify the shape of the underlying hazard. The hazard function \(h(t)\) and survival function \(S(t)\) are thus specified fully by a set of parameters.

3 Why parametric survival models?

Flexibility vs. specification

The choice of a survival model often involves a trade-off between being “data-driven” (flexible) and “model-driven” (specified).

  • Non-parametric methods (e.g., Kaplan-Meier): These make no assumptions about the distribution of survival times. While flexible, they result in “step functions” that are tied strictly to the observed event times.
  • Semi-parametric methods (e.g., Cox proportional hazards model): The Cox model assumes proportional hazards but leaves the shape of the baseline hazard unspecified. This makes it widely applicable and robust, but it cannot provide a full picture of the hazard over time.
  • Parametric models: These assume a specific functional form (shape) for the baseline hazard. While this requires more justification, it offers several powerful advantages.

Key advantages of parametric models

  • Completeness: Both the hazard function \(h(t)\) and the survival function \(S(t)\) are fully specified for all values of \(t\).
  • Smoothness and continuity: Unlike non-parametric step functions (which only change at the exact moment of an observed event), parametric models provide smooth, continuous functions. This is often more consistent with biological theory, where the underlying risk (hazard) of a disease usually evolves continuously over time rather than in discrete jumps.
  • Extrapolation: Perhaps the most significant advantage is the ability to predict survival beyond the last observed event.
  • Efficiency: if the chosen distribution is a good fit for the data, parametric models provide more precise inferences (narrower confidence intervals) than semi-parametric models.

4 Common distributions in parametric survival models

Distributions commonly used for parametric survival models include:

  • Exponential: Assumes a constant hazard over time (e.g., a speedometer stuck at 100 km/h)
  • Weibull: A monotonic hazard that can either increase or decrease over time (e.g., a car constantly accelerating or decelerating)
  • Log-logistic/lognormal: Useful for unimodal non-monotic hazard that first increases and then decreases (e.g., recovery from surgery)
  • Generalized gamma: A flexible distribution that can accommodate various different hazard shapes (constant, increasing, decreasing, unimodal, bathtub)
  • Gompertz: Specifically useful for modeling human mortality where hazard increases exponentially with age
  • Piecewise exponential: Acts as a bridge between parametric and Cox models. Uses a step-function approach to modelling flexible hazards by assuming a constant hazard with pre-defined time intervals

5 Exponential distribution

Suppose that the survival times of \(n\) individuals, \(t_1,t_2,…,t_n\) are assumed to have an exponential distribution, then

\[f(t) = \lambda e^{-\lambda t}\]

\[S(t) = P( > t) = \int_t^{\infty} f(u)du = \int_t^{\infty} \lambda e^{-\lambda u} = e^{-\lambda t}\]

\[h(t) = \frac{f(t)}{S(t)} = \lambda, \text{a constant hazard rate}\]

6 Weibull distribution

The constant hazard rate assumption underlying exponentially distributed survival times can be unrealistic.

The Weibull distribution allows more flexibility in the shape of the hazard function:

\[h(t) = \lambda \gamma t^{\gamma-1}\] \[S(t) = \exp[-H(t)] = \exp[-\int_0^t h(u)du] = \exp(-\lambda t^\gamma)\]

\[f(t) = h(t)S(t) = \lambda \gamma t^{\gamma-1} \exp(-\lambda t^\gamma),\] where \(\lambda\) is the scale parameter and \(\gamma\) is the shape parameter.

The Weibull simplifies to the exponential when \(\gamma = 1\)

The scale parameter determines the stretch or spread of the distribution along the time axis, whereas the shape parameter determines if the hazard is increasing, decreasing, or constant. If you increase the scale parameter while keeping the shape constant, you are essentially stretching the survival time, meaning the event happens later on average. In terms of the speedometer analogy, if the shape is how the car accelerates, the scale is the gear the car is in, determining the overall range of speeds possible.

Weibull hazard, survival and probability density functions for different values of the shape parameter while assuming scale parameter is 1 are shown below:

Show R code
# make sure you have installed the tidyverse package
# install.packages("tidyverse")
# load the tidyverse package
library(tidyverse) # we are using tidyverse tools such as the pipe (%>%) and dplyr

# Weibull hazard, survival and probability density functions

# assumes scale parameter is 1
# four different values for the shape parameter: 0.1, 1, 2, 3
data_weibull <- data.frame(expand.grid(Time = seq(0.01,2,0.01),
                              Gamma = c(0.1, 1, 2, 3))) %>%
  mutate(Hazard = Gamma*Time^(Gamma - 1),
         Survival = exp(-Time^Gamma),
         PDF = Gamma*Time^(Gamma - 1)*exp(-Time^Gamma),
         c = recode_factor(Gamma, `0.1` = "gamma == 0.1", 
                           `1` = "gamma == 1.0", 
                           `2` = "gamma == 2", 
                           `3` = "gamma == 3"))

# Weibull hazard functions
ggplot(data_weibull, aes(x=Time, y=Hazard, group=c)) +
  geom_line(aes(linetype=c, color=c), linewidth = 0.8) +
  labs(color = '', linetype = '') +
  scale_colour_discrete(labels = scales::parse_format()) +
  scale_linetype_discrete(labels = scales::parse_format()) +
  theme_bw()

Show R code
# Weibull survival functions
ggplot(data_weibull, aes(x=Time, y=Survival, group=c)) +
  geom_line(aes(linetype=c, color=c), linewidth = 0.8) +
  labs(color = '', linetype = '') +
  scale_colour_discrete(labels = scales::parse_format()) +
  scale_linetype_discrete(labels = scales::parse_format()) +
  theme_bw()

Show R code
# Weibull probability density functions
ggplot(data_weibull, aes(x=Time, y=PDF, group=c)) +
  geom_line(aes(linetype=c, color=c), linewidth = 0.8) +
  labs(color = '', linetype = '', y="Probability Density Function") +
  scale_colour_discrete(labels = scales::parse_format()) +
  scale_linetype_discrete(labels = scales::parse_format()) +
  theme_bw()

A unique property of the Weibull distribution is that regardless of the shape, the scale parameter represents the time point by which approximately 63.2% of the population is expected to have experienced the event. If the scale = 1, it means that at 1 unit of time (e.g., 1 year or 1 month, depending on your data), roughly 63% of your subjects will have had the event, as can be seen from the survival distribution above.

7 Fitting a parametric survival model

Parametric survival models are typically fitted to observed survival data using the method of maximum likelihood. We construct a likelihood function that represents the probability of observing our specific data given a set of parameters. Because the likelihood equations for survival data are often complex, in practice a numerical procedure - most commonly the Newton-Raphson algorithm - is used to find the values of parameters that maximise the likelihood function.

8 Data example: The Gehan leukaemia data

To demonstrate fitting parametric survival models, we will use a data example by Gehan (1965) (which originated from a randomised controlled trial by Freireich et al. 1963), in which

  • 42 children with acute leukaemia were enrolled in a randomised controlled trial, in which 21 patients received an experimental treatment (6-mercaptopurine or 6-MP) and the other 21 received a placebo
  • Participants had complete or partial remission status at randomisation, and were matched by remission status (one randomised to 6-MP and the other to placebo)
  • The participants were then followed until their leukaemia returned (relapse) or the end of study (administrative censoring)

The Gehan study data is included in the MASS package:

Show R code
# make sure you have installed the MASS package
# install.packages("MASS")
# load the MASS package
library(MASS)

# View gehan data for the leukaemia trial in MASS package gehan
MASS::gehan
   pair time cens   treat
1     1    1    1 control
2     1   10    1    6-MP
3     2   22    1 control
4     2    7    1    6-MP
5     3    3    1 control
6     3   32    0    6-MP
7     4   12    1 control
8     4   23    1    6-MP
9     5    8    1 control
10    5   22    1    6-MP
11    6   17    1 control
12    6    6    1    6-MP
13    7    2    1 control
14    7   16    1    6-MP
15    8   11    1 control
16    8   34    0    6-MP
17    9    8    1 control
18    9   32    0    6-MP
19   10   12    1 control
20   10   25    0    6-MP
21   11    2    1 control
22   11   11    0    6-MP
23   12    5    1 control
24   12   20    0    6-MP
25   13    4    1 control
26   13   19    0    6-MP
27   14   15    1 control
28   14    6    1    6-MP
29   15    8    1 control
30   15   17    0    6-MP
31   16   23    1 control
32   16   35    0    6-MP
33   17    5    1 control
34   17    6    1    6-MP
35   18   11    1 control
36   18   13    1    6-MP
37   19    4    1 control
38   19    9    0    6-MP
39   20    1    1 control
40   20    6    0    6-MP
41   21    8    1 control
42   21   10    0    6-MP

When you load the dataset in R, you will see the following variables in the data:

  • pair: Label for matched pair (we will ignore this aspect in our examples)

  • time: The time (in weeks) in remission until the return of the disease (relapse)

  • cens: The censoring indicator (0 = censoring, 1 = relapse)

  • treat: The treatment arm (control, 6-MP).

9 Modelling approaches

To evaluate how covariates affect survival times, two following modelling approaches can be employed:

1. Proportional hazards (PH) models

A significant benefit of the direct mathematical link between the hazard function and the survival function is that any covariate effect modelled on the hazard scale translates intuitively to the survival scale. Consequently, regression models in survival analysis are typically specified at the level of the hazard.

The most common approach involves relative risk models, which assume that covariates have a multiplicative effect on the hazard. For an individual \(i\) (where \(i = 1, \dots, n\)) with a set of \(p\) covariates \(x_i = (x_{i1}, \dots, x_{ip})\), the proportional hazards (PH) model defines the hazard function as:

\[h_i(t) = h(t | x_i) = h_0(t) \cdot \psi(x_i),\]

where

  • \(h_0(t)\) is the baseline hazard, a time-dependent function representing the risk for an individual with all covariate values set to zero. This baseline is common to all individuals in the study.
  • \(\psi(x_i)\) is the multiplier, a non-negative function of the covariates that scales the baseline hazard up or down.

To ensure the hazard remains positive, we typically use the exponential of a linear combination of the covariates (\(x_i \beta\)). The resulting model is expressed as:

\[h_i(t) = h(t | x_i) = h_0(t) \exp(x_i \beta)\] On the log-hazard scale this model writes down

\[\log h_i(t) = \log h(t | x_i) = \log h_0(t) + x_i \beta,\] showing how the covariates shift the hazard proportionally, as the distance between the log-hazard curves remains constant over time.

In our data example, we would be interested to know “Does the treatment lower the risk of leukaemia relapse compared to the placebo?”

The results from proportional hazards models are typically reported as hazard ratios (HR). The main assumption of this model is the proportional hazards assumption, i.e., that the ratio of the hazards for two subjects with covariate values \(x_i\) and \(x_j\) is constant over time:

\[ \frac{h(t|x_i)}{h(t|x_j)} = \frac{\exp(x_i\beta)}{\exp(x_j\beta)}\]

2. Accelerated failure time (AFT) models

While the proportional hazards (PH) model is the standard in survival analysis, the accelerated failure time (AFT) model serves as a powerful alternative. It is particularly useful when the proportional hazards assumption is violated or when the research question focuses on the timing of an event rather than the risk of its occurrence.

Because survival time \(T\) is a positive variable, the AFT model is naturally framed as a log-linear regression. This allows us to model the impact of covariates directly on the event time \(T\) using a log link function:

\[\log(T_i) = \mu + x_i\alpha + \sigma\epsilon_i,\]

where

  • \(\mu\) and \(\sigma\) represent the location and scale parameters, respectively
  • \(x_i\alpha\) is the linear predictor of your covariates
  • \(\epsilon_i\) is s the error term, typically assumed to follow a specific parametric distribution (such as Extreme Value, Normal, or Logistic)

The intuition behind the name accelerated failure time model is best understood through the survival function. For any individual \(i\), the AFT model assumes that covariates act directly on the time scale itself:

\[S_i(t) = S_0\left( \phi \cdot t \right),\] where \(S_0(.)\) is the baseline survival function and \(\phi >0\) is an acceleration factor (AF).

More generally, one can write \(\phi = exp(-x_i\alpha)\):

\[S_i(t) = S_0\left( \exp(-x_i\alpha) \cdot t \right)\] In this framework, the covariates act as an accelerator factor \((\phi = \exp(-x_i\alpha))\):

  • Accelerating time, if \(\alpha > 1\), making the time pass ‘faster’ for that individual, so that they reach the event sooner than the baseline
  • Decelerating time, if \(\alpha < 1\), making time ‘slow down’ and extending the survival period compared to the baseline

In AFT model, the covariates thus stretch out or shrink survival time. In our data example, we would be interested to know “How much longer can a patient expect to remain in remission on the treatment compared to the placebo?”

In our data example, we would be interested to know “Does the treatment lower the risk of leukaemia relapse compared to the placebo?”

The results from accelerated failure models are typically reported as time ratios (TR). The main assumption of the AFT model is that the effect of a covariate is to multiply the time scale, i.e., that the ratio of event times in two groups, i.e., exposed (E) and control (C)

\[ TR = \frac{T_E}{T_C} \] is constant over time.

AFT models are usually fully parametric.

In the following we are going to show how to fit various fully parametric proportional hazards and accelerated failure time models in our example dataset.

10 Data example: Fitting an Exponential model

For exponential model, applied to the data example, where \(X\) = group (1 = treatment, 0 = placebo):

1. Proportional hazards (PH) model

\[h(t | X) = h_0(t) \exp(X \beta_1) = \lambda \exp(X \beta_1) = \exp(\beta_0 + X \beta_1) \]

2. Accelerated failure time (AFT) model

\[\log(T) = \mu + X\alpha_1 + \sigma\epsilon,\] where for exponential model, the scale parameter \(\sigma\) is fixed at 1.

Exponential PH and AFT models are thus the same model, the only difference is in their parametrisation.

Fitting exponential proportional hazards (PH) model

We can fit the exponential PH model using flexsurvreg function from the flexsurv package, specifying dist="exp":

Show R code
# make sure you have installed the flexsurv package
# install.packages("flexsurv")
# load the flexsurv package
# note that when you load the flexsurv package, it also automatically loads the 
# survival package, as it is built directly on top of the survival package's 
# infrastructure, and needs the Surv() function and the underlying survival object 
# logic to function
library(flexsurv) # for flexsurvreg for parametric models
library(MASS) # reloaded to be able to run this code block independently

# read in the gehan data for the leukaemia trial in MASS package gehan
leukdata <- MASS::gehan

# Setting control group as the reference group
leukdata$treat <- factor(leukdata$treat, levels = c("control", "6-MP"))

# EXPONENTIAL PH MODEL

# Fitting an exponential PH model using flexsurv package
exp_PH <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "exp")
# Summary shows results in terms of rate by default
# gives rate for the reference group (placebo)
# gives HR
exp_PH
Call:
flexsurvreg(formula = Surv(time, cens) ~ treat, data = leukdata, 
    dist = "exp")

Estimates: 
           data mean  est      L95%     U95%     se       exp(est)  L95%   
rate            NA     0.1154   0.0752   0.1770   0.0252       NA        NA
treat6-MP   0.5000    -1.5266  -2.3075  -0.7457   0.3984   0.2173    0.0995
           U95%   
rate            NA
treat6-MP   0.4744

N = 42,  Events: 30,  Censored: 12
Total time at risk: 541
Log-likelihood = -108.524, df = 2
AIC = 221.0481
Show R code
# This will calculate the actual hazard rate for each group
summary(exp_PH, newdata = data.frame(treat = c("control", "6-MP")), type = "hazard")
treat=control 
   time       est        lcl       ucl
1     1 0.1153846 0.07563433 0.1712061
2     2 0.1153846 0.07563433 0.1712061
3     3 0.1153846 0.07563433 0.1712061
4     4 0.1153846 0.07563433 0.1712061
5     5 0.1153846 0.07563433 0.1712061
6     6 0.1153846 0.07563433 0.1712061
7     7 0.1153846 0.07563433 0.1712061
8     8 0.1153846 0.07563433 0.1712061
9     9 0.1153846 0.07563433 0.1712061
10   10 0.1153846 0.07563433 0.1712061
11   11 0.1153846 0.07563433 0.1712061
12   12 0.1153846 0.07563433 0.1712061
13   13 0.1153846 0.07563433 0.1712061
14   15 0.1153846 0.07563433 0.1712061
15   16 0.1153846 0.07563433 0.1712061
16   17 0.1153846 0.07563433 0.1712061
17   19 0.1153846 0.07563433 0.1712061
18   20 0.1153846 0.07563433 0.1712061
19   22 0.1153846 0.07563433 0.1712061
20   23 0.1153846 0.07563433 0.1712061
21   25 0.1153846 0.07563433 0.1712061
22   32 0.1153846 0.07563433 0.1712061
23   34 0.1153846 0.07563433 0.1712061
24   35 0.1153846 0.07563433 0.1712061

treat=6-MP 
   time        est        lcl        ucl
1     1 0.02506964 0.01277963 0.04621725
2     2 0.02506964 0.01277963 0.04621725
3     3 0.02506964 0.01277963 0.04621725
4     4 0.02506964 0.01277963 0.04621725
5     5 0.02506964 0.01277963 0.04621725
6     6 0.02506964 0.01277963 0.04621725
7     7 0.02506964 0.01277963 0.04621725
8     8 0.02506964 0.01277963 0.04621725
9     9 0.02506964 0.01277963 0.04621725
10   10 0.02506964 0.01277963 0.04621725
11   11 0.02506964 0.01277963 0.04621725
12   12 0.02506964 0.01277963 0.04621725
13   13 0.02506964 0.01277963 0.04621725
14   15 0.02506964 0.01277963 0.04621725
15   16 0.02506964 0.01277963 0.04621725
16   17 0.02506964 0.01277963 0.04621725
17   19 0.02506964 0.01277963 0.04621725
18   20 0.02506964 0.01277963 0.04621725
19   22 0.02506964 0.01277963 0.04621725
20   23 0.02506964 0.01277963 0.04621725
21   25 0.02506964 0.01277963 0.04621725
22   32 0.02506964 0.01277963 0.04621725
23   34 0.02506964 0.01277963 0.04621725
24   35 0.02506964 0.01277963 0.04621725

The output for the intercept is \(\exp(\hat{\beta_0})\), giving the baseline hazard for the reference group (placebo) directly: \(\hat{\lambda_0}=0.115\). \(\hat{\beta_1} = -1.527\), on the other hand, is reported on the log-hazard scale, and the sign of the raw coefficient tells you the direction of the effect immediately (e.g., negative is a reduction in risk). However, the output also gives you the \(\exp(\hat{\beta_1})\), which is the hazard ratio (HR) we are typically interested in:

\[HR = \frac{\exp(\hat{\beta_0}+\hat{\beta_1})}{\exp(\hat{\beta_0})} = \exp(\hat{\beta_1}) = \exp(-1.527) = 0.217\] We can also calculate the hazard for the treatment group \(\hat{\lambda_1} = \exp(\hat{\beta_0}+\hat{\beta_1}) = 0.115*0.217 = 0.025\)

Assuming the data follow exponential distribution, the treatment lowers the risk of leukaemia relapse by 78% compared to the placebo (HR = 0.217, 95% CI = 0.10, 0.474).

Fitting exponential accelerated failure time (AFT) model

We can fit the exponential accelerated failure time model using survreg function from the survival package, specifying dist="exponential":

Show R code
# make sure you have installed the survival package
# install.packages("survival")
# load the survival package
# note that when you loaded the flexsurv package above, it also automatically 
# loaded the survival package, as it is built directly on top of the survival 
# package's infrastructure, and needs the Surv() function and the underlying 
# survival object logic to function
library(survival) # for survreg for parametric models

# EXPONENTIAL AFT MODEL

# Fitting an exponential AFT model using survreg package
exp_AFT <- survreg(Surv(time,cens) ~ treat, data = leukdata, dist = "exponential")
summary(exp_AFT)

Call:
survreg(formula = Surv(time, cens) ~ treat, data = leukdata, 
    dist = "exponential")
            Value Std. Error    z       p
(Intercept) 2.159      0.218 9.90 < 2e-16
treat6-MP   1.527      0.398 3.83 0.00013

Scale fixed at 1 

Exponential distribution
Loglik(model)= -108.5   Loglik(intercept only)= -116.8
    Chisq= 16.49 on 1 degrees of freedom, p= 4.9e-05 
Number of Newton-Raphson Iterations: 4 
n= 42 
Show R code
# Extract coefficients
coeffs <- coef(exp_AFT)
coeffs
(Intercept)   treat6-MP 
   2.159484    1.526614 
Show R code
# Mean survival time in the control (placebo) group
T_C <- exp(coeffs[1])
T_C
(Intercept) 
   8.666667 
Show R code
# Mean survival time in the exposed (treatment) group
T_E <- exp(coeffs[1]+coeffs[2])
T_E
(Intercept) 
   39.88889 
Show R code
# Calculate Time Ratio (TR) for the treatment
time_ratio <- exp(coeffs[2])
time_ratio
treat6-MP 
 4.602564 
Show R code
# Calculate 95% confidence intervals for the TR
conf_intervals <- exp(confint(exp_AFT)[2, ])
conf_intervals
    2.5 %    97.5 % 
 2.108012 10.049088 

As the AFT model is parameterised as a linear regression of log-survival time, the intercept (\(\mu\)) represents the log of the mean survival time for the control (placebo) group, with mean survival time thus \(T_C = \exp(\hat{\mu}) = 8.7\) weeks. The coefficient \(\hat{\alpha_1} = 1.527\) is also expressed on the log-time scale, and the sign of the raw coefficient tells you the direction of the effect immediately (e.g., positive means the treatment increases the time in remission). We can calculate the mean survival time for the exposed (treatment) group \(T_E = \exp(\hat{\mu}+\hat{\alpha_1}) = \exp(2.159+1.527) = 39.9\) weeks. We can also calculate the time ratio (TR) that we are typically interested in:

\[TR = \frac{T_E}{T_C} = \frac{\exp(\hat{\mu}+\hat{\alpha_1})}{\hat{\exp(\hat{\mu})}} = \exp(\hat{\alpha_1}) = \exp(1.527) = 4.6\] Thus, assuming the data follow exponential distribution, the treatment extends the time in remission 4.6 times (TR = 4.6, 95% CI = 2.1, 10.0).

Note that for exponential model the coefficient \(\alpha_1\) for the AFT model is the opposite to the coefficient \(\beta_1\) in PH model: \(\beta_1 = - \alpha_1\). The hazard ratio for exponential PH model can thus also be obtained from the survreg output as \(\exp(-1.527) = 0.22\).

We can also use the predict() function from the stats package (which is loaded by default) to compute for example median survival time . When a survreg object is used, a specific predict.survreg() function is called:

Show R code
# EXPONENTIAL AFT MODEL

# Function predict can compute median survival time 
# Median survival time for placebo and treatment group
median_times <- predict(exp_AFT, newdata = data.frame(treat = c("control", "6-MP")), 
                        type = "quantile", p = 0.5)
median_times
        1         2 
 6.007276 27.648871 
Show R code
# predict() function does not work for flexsurv objects but we can use summary
median_by_group <- summary(exp_PH, type = "quantile", quantiles = 0.5)
# stack the list elements into one table
do.call(rbind, median_by_group)
              quantile       est       lcl       ucl
treat=control      0.5  6.007276  3.893135  9.297986
treat=6-MP         0.5 27.648871 14.355990 54.687991
Show R code
# similarly for mean
mean_by_group <- summary(exp_PH, type = "mean")
do.call(rbind, mean_by_group)
                    est       lcl      ucl
treat=control  8.666667  5.804312 13.69279
treat=6-MP    39.888889 20.725716 76.59898
Show R code
# can also be used to obtain restricted mean survival time to chosen tau = t
#summary(exp_PH, type = "rmst", t = 10)

Note that the ratio of the median survival times \(TR = \frac{27.6}{6.0} = 4.6\).

11 Data example: Fitting a Weibull model

For Weibull model, applied to the data example, where \(X\) = group (1 = treatment, 0 = placebo):

1. Proportional hazards (PH) model

\[h(t | X) = h_0(t) \exp(X \beta_1) = \lambda \gamma t^{\gamma-1} \exp(X \beta_1) = \exp(\beta_0 + X \beta_1)\gamma t^{\gamma-1}, \]

where \(h_0(t) = \exp(\beta_0)\gamma t^{\gamma-1}\)

2. Accelerated failure time (AFT) model

\[\log(T) = \mu + X\alpha_1 + \sigma\epsilon,\] The following relationships apply:

\(\lambda = \exp(\frac{-\mu}{\sigma}), \gamma = \frac{1}{\sigma}, \beta_1 = \frac{-\alpha_1}{\sigma}\)

For further detail and proof see:

Catherine Legrand. Advanced survival models. 2021.

Note also that the Weibull distribution is mathematically unique in that it is the only distribution that satisfies both the proportional hazards (PH) and accelerated failure time (AFT) assumptions simultaneously.

Fitting Weibull proportional hazards (PH) model

We can fit the Weibull PH model using flexsurvreg function from the flexsurv package, specifying dist="weibullPH":

Show R code
# WEIBULL PH MODEL

# Fitting a Weibull PH model using flexsurv package
wei_PH <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "weibullPH")

# gives estimate for shape parameter and scale parameter
# gives HR
wei_PH
Call:
flexsurvreg(formula = Surv(time, cens) ~ treat, data = leukdata, 
    dist = "weibullPH")

Estimates: 
           data mean  est      L95%     U95%     se       exp(est)  L95%   
shape           NA     1.3658   1.0233   1.8228   0.2012       NA        NA
scale           NA     0.0464   0.0155   0.1385   0.0259       NA        NA
treat6-MP   0.5000    -1.7309  -2.5405  -0.9212   0.4131   0.1771    0.0788
           U95%   
shape           NA
scale           NA
treat6-MP   0.3980

N = 42,  Events: 30,  Censored: 12
Total time at risk: 541
Log-likelihood = -106.5795, df = 3
AIC = 219.159
Show R code
# This will calculate the actual hazard rates for each group at different times
summary(wei_PH, newdata = data.frame(treat = c("control", "6-MP")), type = "hazard")
treat=control 
   time        est        lcl       ucl
1     1 0.06335542 0.02896587 0.1482740
2     2 0.08163713 0.04720033 0.1581785
3     3 0.09468789 0.06049797 0.1641259
4     4 0.10519418 0.06981734 0.1734936
5     5 0.11413984 0.07694147 0.1813703
6     6 0.12201084 0.08232592 0.1947539
7     7 0.12908765 0.08620431 0.2094615
8     8 0.13554880 0.08851323 0.2256332
9     9 0.14151586 0.09128896 0.2474927
10   10 0.14707581 0.09341675 0.2647655
11   11 0.15229335 0.09504637 0.2864100
12   12 0.15721804 0.09713373 0.3022444
13   13 0.16188883 0.09785789 0.3142957
14   15 0.17058780 0.10088346 0.3459331
15   16 0.17466250 0.10215017 0.3646424
16   17 0.17857871 0.10300762 0.3837658
17   19 0.18599338 0.10442243 0.4148969
18   20 0.18951572 0.10481474 0.4312832
19   22 0.19623882 0.10644751 0.4632220
20   23 0.19945546 0.10686611 0.4778181
21   25 0.20563205 0.10744168 0.5050704
22   32 0.22506277 0.10938770 0.6014698
23   34 0.23010903 0.10977422 0.6263979
24   35 0.23256172 0.10935084 0.6386798

treat=6-MP 
   time        est         lcl        ucl
1     1 0.01122214 0.003491616 0.03682991
2     2 0.01446038 0.005763104 0.03873474
3     3 0.01677206 0.007450192 0.04018123
4     4 0.01863304 0.008742910 0.04164122
5     5 0.02021758 0.010023985 0.04333368
6     6 0.02161177 0.010865688 0.04489062
7     7 0.02286529 0.011974406 0.04622443
8     8 0.02400975 0.012744798 0.04760341
9     9 0.02506669 0.013640158 0.04926894
10   10 0.02605153 0.014455278 0.05044640
11   11 0.02697571 0.014903805 0.05236010
12   12 0.02784802 0.015284642 0.05474684
13   13 0.02867536 0.015707240 0.05696528
14   15 0.03021621 0.016631869 0.05993230
15   16 0.03093796 0.016956210 0.06185915
16   17 0.03163164 0.017076900 0.06384246
17   19 0.03294499 0.017426324 0.06844853
18   20 0.03356891 0.017568077 0.07004709
19   22 0.03475977 0.017821656 0.07361532
20   23 0.03532953 0.017919608 0.07563711
21   25 0.03642359 0.018168720 0.08038690
22   32 0.03986535 0.018915943 0.09371726
23   34 0.04075920 0.019138670 0.09827424
24   35 0.04119364 0.019230892 0.10042881

The output gives estimates for the shape \((\hat{\gamma})\) and scale \((\hat{\sigma})\) parameters on the natural scale. The estimate for \((\hat{\gamma})\) is 1.366 suggesting an increase in the hazard as survival time increases (because \(\gamma>1\)). The 95% CI for \(\hat{\gamma}\) (1.023, 1.822) does not include 1, suggesting that the exponential model is not appropriate.

\(\hat{\beta_1} = -1.731\), on the other hand, is reported on the log-hazard scale, and the sign of the raw coefficient tells you the direction of the effect immediately (e.g., negative is a reduction in risk). However, the output also gives you the \(\exp(\hat{\beta_1})\), which is the hazard ratio (HR) we are typically interested in:

\[HR = \frac {\exp(\hat{\beta_0}+\hat{\beta_1}) \hat{\gamma}t^{\hat{\gamma}-1}} {\exp(\hat{\beta_0}) \hat{\gamma}t^{\hat{\gamma}-1}} = \exp(\hat{\beta_1}) = \exp(-1.731) = 0.177\]

Assuming the data follow Weibull distribution, the treatment lowers the risk of leukaemia relapse by 82% compared to the placebo (HR = 0.177, 95% CI = 0.079, 0.398).

Fitting Weibull accelerated failure time (AFT) model

We can fit the Weibull accelerated failure time model using survreg function from the survival package, specifying dist="weibull":

Show R code
# WEIBULL AFT MODEL

# Fitting a Weibull model using survreg package
wei_AFT_surv <- survreg(Surv(time,cens) ~ treat, data = leukdata, dist = "weibull")
summary(wei_AFT_surv)

Call:
survreg(formula = Surv(time, cens) ~ treat, data = leukdata, 
    dist = "weibull")
             Value Std. Error     z       p
(Intercept)  2.248      0.166 13.55 < 2e-16
treat6-MP    1.267      0.311  4.08 4.5e-05
Log(scale)  -0.312      0.147 -2.12   0.034

Scale= 0.732 

Weibull distribution
Loglik(model)= -106.6   Loglik(intercept only)= -116.4
    Chisq= 19.65 on 1 degrees of freedom, p= 9.3e-06 
Number of Newton-Raphson Iterations: 5 
n= 42 
Show R code
# Extract coefficients
coeffs <- coef(wei_AFT_surv)
coeffs
(Intercept)   treat6-MP 
   2.248352    1.267335 
Show R code
# Calculate Time Ratio (TR) for the treatment
time_ratio <- exp(coeffs["treat6-MP"])
time_ratio
treat6-MP 
 3.551374 
Show R code
# Calculate 95% confidence intervals for the TR
conf_intervals <- exp(confint(wei_AFT_surv)["treat6-MP", ])
conf_intervals
   2.5 %   97.5 % 
1.931876 6.528503 
Show R code
# Function predict can compute median survival time 
# Median survival time for placebo and treatment group
median_times <- predict(wei_AFT_surv, newdata = data.frame(treat = c("control", "6-MP")), 
                        type = "quantile", p = 0.5)
median_times
        1         2 
 7.242697 25.721526 

Consistent with the log-linear structure of the AFT model, the intercept/location \((\hat{\mu} = 2.248)\), coefficient \((\hat{\alpha_1} = 1.267)\), and scale \((\hat{\sigma} = -0.312)\) are all reported on the log-time scale, although the actual scale \((\exp(\hat{\sigma}) = 0.732\)) is also given.

From the output for scale, we see that it’s not equal to 1 (or log(Scale) equal to 0), with p-value 0.03. This gives us evidence that the Weibull model fits better than the exponential model (as \(\sigma \ne 1)\) as in exponential model).

The positive sign of the raw coefficient tells us that the treatment increases the time in remission. As before, we can calculate the time ratio (TR) that we are typically interested in:

\[TR = \frac{T_E}{T_C} = \exp(\hat{\alpha_1}) = \exp(1.267) = 3.55\] We can extract 95% confidence intervals for the TR as shown above, and conclude that - assuming the data follow Weibull distribution - the treatment extends the time in remission 3.55 times (TR = 3.55, 95% CI = 1.93, 6.53).

Again, we can also use the predict() function to compute for example median survival times in the two groups, and note that the ratio of the median survival times \(TR = \frac{25.7}{7.24} = 3.55\).

Note that we can also obtain the estimated values in the Weibull PH model parameterisation from the survreg output by following the equations shown before:

\(\lambda = \exp(\frac{-\mu}{\sigma}), \gamma = \frac{1}{\sigma}, \beta_1 = \frac{-\alpha_1}{\sigma}\)

Show R code
lambda <- exp( -wei_AFT_surv$coef[1] / wei_AFT_surv$scale )
lambda
(Intercept) 
 0.04638848 
Show R code
gamma <- 1 / wei_AFT_surv$scale
gamma
[1] 1.365758
Show R code
beta <- -wei_AFT_surv$coef[2] / wei_AFT_surv$scale
beta
treat6-MP 
-1.730872 
Show R code
HR <- exp(beta)
HR
treat6-MP 
0.1771299 

However, there is also a package parfm that does these transformations and allows you to fit the Weibull PH model:

Show R code
# make sure you have installed the parfm package
# install.packages("parfm")
# load the parfm package
library(parfm) 

# WEIBULL PH MODEL

# Fitting a Weibull model using survreg package
wei_AFT_surv2 <- parfm(Surv(time,cens) ~ treat, data = leukdata, dist = "weibull", frailty = "none")
wei_AFT_surv2

Frailty distribution: none 
Baseline hazard distribution: Weibull 
Loglikelihood: -106.579 

          ESTIMATE SE    p-val    
rho        1.366   0.200          
lambda     0.046   0.026          
treat6-MP -1.731   0.413 <.001 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Show R code
summary(wei_AFT_surv2)
    ESTIMATE              SE              p-val         
 Min.   :-1.73087   Min.   :0.02578   Min.   :2.79e-05  
 1st Qu.:-0.84224   1st Qu.:0.11310   1st Qu.:2.79e-05  
 Median : 0.04639   Median :0.20043   Median :2.79e-05  
 Mean   :-0.10624   Mean   :0.21309   Mean   :2.79e-05  
 3rd Qu.: 0.70607   3rd Qu.:0.30675   3rd Qu.:2.79e-05  
 Max.   : 1.36576   Max.   :0.41308   Max.   :2.79e-05  
                                      NA's   :2         
Show R code
HR <- exp(wei_AFT_surv2[c(3)])
HR
[1] 0.1771299
Show R code
HR_CI <- ci.parfm(wei_AFT_surv2, level = 0.05)
HR_CI
            low    up
treat6-MP 0.079 0.398

Finally, we can also fit the Weibull accelerated failure time model using flexsurvreg function from the flexsurv package, specifying dist="weibull":

Show R code
# WEIBULL AFT MODEL

# Fitting a Weibull AFT model using flexsurvreg package
wei_AFT <- flexsurvreg( Surv(time,cens) ~ treat, data = leukdata, dist = "weibull")
wei_AFT
Call:
flexsurvreg(formula = Surv(time, cens) ~ treat, data = leukdata, 
    dist = "weibull")

Estimates: 
           data mean  est     L95%    U95%    se      exp(est)  L95%    U95%  
shape          NA      1.366   1.023   1.823   0.201      NA        NA      NA
scale          NA      9.472   6.842  13.114   1.572      NA        NA      NA
treat6-MP   0.500      1.267   0.658   1.876   0.311   3.551     1.932   6.529

N = 42,  Events: 30,  Censored: 12
Total time at risk: 541
Log-likelihood = -106.5795, df = 3
AIC = 219.159

The output for the shape parameter (\(\hat{\gamma}\)) is the same as in the Weibull PH model (reported on the natural scale). As mentioned before, a unique property of the Weibull distribution is that regardless of the shape, the scale parameter represents the time point by which approximately 63.2% of the population is expected to have experienced the event. In the AFT parameterisation of the Weibull model, the scale parameter in the output refers to this time point, i.e. to time \(t\) at which \(S(t)≈0.368\), which in this data example is 9.472 weeks for relapse in the reference (placebo) group.

Consistent with the log-linear structure of the AFT model, the covariate effect \(\hat{\alpha_1} = 1.267\) is reported on the log-time scale. Similarly as shown above for the survreg output, we can calculate the time ratio (TR) that we are typically interested in:

\[TR = \frac{T_E}{T_C} = \exp(\hat{\alpha_1}) = \exp(1.267) = 3.55\] Note, however, that the flexsurvreg output also automatically provides the TR and its 95% confidence interval.

12 Model validation

Before interpreting the coefficients or baseline parameters from a parametric survival model, we must verify that the chosen distribution is a reasonable approximation of the underlying data-generating process. Unlike semi-parametric models, parametric estimation relies heavily on the validity of the functional form and the behavior of the covariates.

First, we must ensure the distributional form (e.g., Exponential or Weibull) matches the shape of the observed survival history. Second, we must verify the structural assumptions — specifically the proportional hazards (PH) or accelerated failure time (AFT) assumptions. We check this by ensuring that the effect of our covariates remains constant across the entire follow-up period.

1. Visual comparison with non-parametric estimates

The simplest check is to plot the fitted parametric survival curve \(\hat{S}(t)\) directly over the Kaplan-Meier (KM) curve, that is ‘data-driven’ and makes no functional form assumptions. If the smooth parametric curve deviates significantly from the KM steps, the model is likely misspecified.

KM versus Exponential fit:

Show R code
# 1. Fit the Kaplan-Meier (observed) by treatment group
km <- survfit(Surv(time, cens) ~ treat, data=leukdata)

# 2. Fit the parametric exponential survival model (predicted) by treatment group
exp_AFT <- survreg(Surv(time,cens) ~ treat, data = leukdata, dist = "exponential")

# 3. Visual suitability check: survival by treatment group

# Plot the KM curves
plot(km, col = c("red", "blue"), conf.int = FALSE, 
     main = "KM vs Exponential fit", xlab = "Weeks", ylab = "Survival Probability")
legend("bottomleft", legend = levels(leukdata$treat), col = c("red", "blue"), lty = 1:2)

# Add the Exponential theoretical curves
times <- seq(0, max(leukdata$time), length.out = 100)

# Calculate S(t) for Group 1 (reference = placebo)
# lambda = exp(-intercept)
lambda_ref <- exp(-coef(exp_AFT)[1])
lines(times, exp(-lambda_ref * times), col = "red", lwd = 2, lty=2)

# Calculate S(t) for Group 2 (treatment)
# lambda = exp(-(intercept + coefficient)
lambda_treat <- exp(-(exp_AFT$coef[1] + exp_AFT$coef[2]))
lines(times, exp(-lambda_treat * times), col = "blue", lwd = 2, lty = 2)

KM versus Exponential and Weibull fit:

Show R code
# 1. Fit the Kaplan-Meier (observed) by treatment group
# km <- survfit(Surv(time, cens) ~ treat, data=leukdata)

# 2. Fit the parametric exponential survival model (predicted) by treatment group
# exp_AFT <- survreg(Surv(time,cens) ~ treat, data = leukdata, dist = "exponential")

# 3. Fit the parametric Weibull survival model (predicted) by treatment group
wei_AFT_surv <- survreg(Surv(time,cens) ~ treat, data = leukdata, dist = "weibull")

# 4. Visual suitability check: survival by treatment group

# Plot the KM curves
plot(km, col = c("red", "blue"), conf.int = FALSE, 
     main = "KM vs Exponential and Weibull fit", xlab = "Weeks", ylab = "Survival Probability")
legend("bottomleft", legend = levels(leukdata$treat), col = c("red", "blue"), lty = 1:2)

# Add the Exponential and Weibull theoretical curves
times <- seq(0, max(leukdata$time), length.out = 100)

# Exponential
# Calculate S(t) for Placebo (reference)
# lambda = exp(-intercept)
lambda_ref_exp <- exp(-exp_AFT$coef[1])
lines(times, exp(-lambda_ref_exp * times), col = "red", lwd = 2, lty=2)

# Calculate S(t) for 6-MP
# lambda = exp(-(intercept + coefficient_treatment)
lambda_treat_exp <- exp(-(exp_AFT$coef[1] + exp_AFT$coef[2]))
lines(times, exp(-lambda_treat_exp * times), col = "blue", lwd = 2, lty = 2)

# Weibull
# we know S(t) for Weibull model is exp(-lambda*t^gamma)
# Calculate S(t) for Placebo (reference)
# lambda = exp(-intercept/scale)
lambda_ref_wei <- exp(-wei_AFT_surv$coef[1]/wei_AFT_surv$scale)
gamma_wei <- 1/wei_AFT_surv$scale
lines(times, exp(-lambda_ref_wei * times^gamma_wei), col = "red", lwd = 2, lty=3)

# Calculate S(t) for 6-MP
# lambda = exp(-(intercept + coefficient_treatment)/scale
lambda_treat_wei <- exp(-(wei_AFT_surv$coef[1] + wei_AFT_surv$coef[2])/wei_AFT_surv$scale)
lines(times, exp(-lambda_treat_wei * times^gamma_wei), col = "blue", lwd = 2, lty = 3)

We are essentially checking if the parametric curve can adequately capture the steps of the empirical Kaplan-Meier curve. However, the survival function is naturally curved and so visual inspection is not always a good check for model fit.

We could also plot the hazard functions:

Show R code
plot(exp_PH, type = "hazard", main = "Exponential hazard fit")

Show R code
plot(wei_PH, type = "hazard", main = "Weibull hazard fit")

2. Linear transformation (graphical checks)

Alternatively, we can exploit the mathematical properties of these distributions by transforming the survival function into a linear form. If the distribution fits, the transformed data should fall on a straight line.

For exponential model, the survival function is \(S(t)=\exp(-\lambda t)\) and \(-\log{S(t)}=-\lambda t\), and thus the cumulative hazard, \(-\log{S(t)}\), is linearly related over time with \(\lambda\).

We can thus check whether plotting the cumulative hazard against time gives a straight line starting at the origin:

Show R code
# make sure you have installed the survminer package
# install.packages("survminer")
# load the survminer package
library(survminer) # for surv_summary() function
library(tidyverse) # reloaded to be able to run this block independently
library(survival) # reloaded to be able to run this block independently

# Fitting cumulative hazard against time

# Kaplan-Meier (observed) by treatment group
km <- survfit(Surv(time, cens) ~ treat, data=leukdata)
cum_data <- km %>% surv_summary(data = leukdata)

cum_data <- cum_data %>%
  mutate(ch = -log(surv))

# using coefficients from this model:
# exp_AFT <- survreg( Surv(time,cens) ~ treat, data = leukdata, dist = "exponential")
ggplot(cum_data, aes(x = time, y = ch, color = treat)) +
  geom_line() +
  geom_point() +
  # we know -logS(t) = lambda * t
  # intercept: 0
  # slope: lambda 
  geom_abline(aes(intercept = 0,
                  slope = exp(-exp_AFT$coef[1])), linetype=2) +
  geom_abline(aes(intercept = 0,
                  slope = exp(-(exp_AFT$coef[1]+exp_AFT$coef[2]))), linetype=2) +
  theme_bw()

If the hazard is constant, as the exponential model assumes, the plot should be a straight diagonal line. Both plots look reasonably straight suggesting that the exponential assumption is reasonable.

We can also use this plot to assess proportional hazards assumption in the exponential PH model: If the hazards are proportional, the lines will diverge at a constant rate. Proportionality is thus visualised as two straight lines that do not cross and stay the same distance apart. That appears to be true here.

For Weibull model, the survival function is \(S(t)=\exp(-\lambda t^\gamma)\) and \(\log(-\log{S(t)}) = \log \lambda + \log t\), and thus the log cumulative hazard, \(\log(-\log{S(t)})\), is linear with log of time.

We can thus check whether plotting the log cumulative hazard against log of time gives a straight line:

Show R code
# Fitting log cumulative hazard against time

# Kaplan-Meier (observed) by treatment group
km <- survfit( Surv(time, cens) ~ treat, data=leukdata)
cum_data <- km %>% surv_summary(data = leukdata) 

cum_data <- cum_data %>%
  mutate(log_ch = log(-log(surv)))

# using coefficients from this model:
# wei_AFT_surv <- survreg( Surv(time,cens) ~ treat, data = leukdata, dist = "weibull")
ggplot(cum_data, aes(x = log(time), y = log_ch, color = treat)) +
  geom_line() +
  geom_point() +
  # we know log(-logS(t) = log(lambda)+gamma*log(t)
  # intercept: log(lambda) = log * exp(-coef/scale)
  # slope: gamma = 1/scale
  geom_abline(aes(intercept = -wei_AFT_surv$coef[1]/wei_AFT_surv$scale,
                  slope = 1/wei_AFT_surv$scale), linetype=2) +
  geom_abline(aes(intercept = -(wei_AFT_surv$coef[1]+wei_AFT_surv$coef[2])/wei_AFT_surv$scale,
                  slope = 1/wei_AFT_surv$scale), linetype=2) +
  theme_bw()

Both plots look reasonably straight suggesting that the Weibull assumption is reasonable.

We can also use this plot to assess the proportional hazards assumption in the Weibull PH model: If the hazards are proportional, the lines will be parallel, i.e., have the same slope (\(\gamma\)). Here, the lines appear parallel.

Note also that for Weibull model, a constant vertical shift (multiplying the hazard) is mathematically indistinguishable from a constant horizontal shift (multiplying the survival time). That is, the Weibull distribution satisfies both the PH and AFT assumptions simultaneously. Thus, we can also conclude that the *accelerated failure time assumption** also holds.

From the log-log plot, we can actually deduct the following things:

  • 1. Parallel straight lines: Weibull and PH (and AFT) assumptions hold
  • 2. Parallel straight lines with slope 1: Exponential model is a good fit
  • 3. Parallel but not straight lines: PH but not Weibull (could use the Cox model)
  • 4. Not parallel and not straight: Weibull is not a good fit, PH violated
  • 5. Not parallel but straight: Weibull is a good fit, but PH and AFT violated, indicates different \(\gamma\) in the two groups

3. Graphical check for accelerated failure time assumption

The accelerated failure time (AFT) model assumes that a covariate (e.g., treatment) multiplies survival time by a constant factor.

We can check whether there is a constant horizontal shift between survival curves plotted against log of time:

Show R code
# make sure you have installed the broom package
# install.packages("broom")
# load the broom package
library(broom) # for tidy() function
library(tidyverse) # reloaded to be able to run this block independently
library(flexsurv) # reloaded to be able to run this block independently

# Plotting survival curves on log time scale to see the AFT model constant shift

# 1. Fit the Kaplan-Meier (observed) by treatment group
km <- survfit( Surv(time, cens) ~ treat, data=leukdata)
#summary(km)
#names(km)

# 2. Create a data frame for ggplot
plot_data <- broom::tidy(km)
#plot_data
# survival probability is now named estimate
#names(plot_data)

# 3. fit Exponential and Weibull
exp_PH <- flexsurvreg( Surv(time, cens) ~ treat, data = leukdata, dist = "exp")
wei_AFT <- flexsurvreg( Surv(time,cens) ~ treat, data = leukdata, dist = "weibull")

# 4. Generate predicted survival probabilities for a range of times
# This creates a smooth curve for the log-time axis
time_seq <- seq(0.1, 35, length.out = 100)
pred_exp <- summary(exp_PH, t = time_seq, tidy = TRUE)
pred_wei <- summary(wei_AFT, t = time_seq, tidy = TRUE)

# 5. Plot everything with a Log-transformed X-axis
ggplot() +
  # Kaplan-Meier steps (the raw data)
  geom_step(data = plot_data, aes(x = time, y = estimate, color = strata), linewidth = 0.8) +
   # Exponential fitted curves
  geom_line(data = pred_exp, aes(x = time, y = est, color = treat), linetype = "dashed", linewidth = 1) +
  # Weibull fitted curves
  geom_line(data = pred_wei, aes(x = time, y = est, color = treat), linetype = "dotted", linewidth = 1) +
  scale_x_log10() + 
  labs(title = "KM vs. Parametric fits on log time scale",
       subtitle = "Dashed = Exponential, Dotted = Weibull",
       x = "Log(Time)", y = "Survival Probability") +
  theme_bw()

The shift between survival curves does seem constant both for the exponential and weibull AFT model.

13 Model selection

Exponential vs Weibull

1. Likelihood Ratio Test (LRT)

Because the Exponential distribution is a special case of the Weibull distribution (where the shape parameter \(\gamma = 1\)), we can use a formal Likelihood Ratio Test (LRT) to determine which model provides a better fit for the data.

To test whether a Weibull \((\lambda, \gamma)\) distribution gives a significantly better fit than an Exponential \((\lambda\)) distribution, we test the following null hypothesis:

\[H_0: \gamma = 1 \quad \text{(Exponential model)}\] \[H_1: \gamma \neq 1 \quad \text{(Weibull model)}\] Under the null hypothesis \(H_0\), the likelihood ratio test statistic \(W\) follows a chi-squared distribution with 1 degree of freedom (\(\chi^2_1\)):

\[W(\hat{\lambda}, 1) = 2 \left[ \log L(\hat{\lambda}, \hat{\gamma}) - \log L(\hat{\lambda_0}, 1) \right] \sim \chi_1^2,\] where \(\hat{\lambda}\) and \(\hat{\gamma}\) are the maximum likelihood estimates (MLEs) for a Weibull model with unconstrained \(\hat{\gamma}\), and \(\hat{\lambda}_0\) is the MLE for a Weibull model with the constraint \(\gamma = 1\) (i.e., MLE for the exponential model).

The Exponential and Weibull PH and AFT model outputs we ran before all provide the necessary log-likelihood values for carrying out the LRT. We can thus carry out the LRT for example in the following way:

Show R code
# Calculate manually

# Using Survreg output above that provides the log-likehood values
# The necessary values can be extracted in the following way:
# the first value of loglik element gives the log-likelihood for the model with only an intercept
# the second value gives the log-likelihood for the model with the covariates
logL_exp <- exp_AFT$loglik[2]
logL_exp
[1] -108.524
Show R code
logL_wei <- wei_AFT_surv$loglik[2]
logL_wei
[1] -106.5795
Show R code
# Calculate the test statistic
LRT_statistic <- 2 * (logL_wei - logL_exp)
LRT_statistic
[1] 3.889116
Show R code
# Calculate the P-value
# We use the Chi-squared distribution with 1 degree of freedom
p_value <- pchisq(LRT_statistic, df = 1, lower.tail = FALSE)
p_value
[1] 0.0486

A significant p-value (\(p < 0.05\)) indicates that allowing the shape parameter \(\gamma\) to vary (Weibull model) significantly improves the model’s log-likelihood compared to fixing it at 1 (Exponential model). We would thus reject the null hypothesis at 5% significance level, and conclude that a general Weibull distribution is better than an exponential distribution in this case.

We could also use the anova() function from the stats package to compare two nested survreg objects:

Show R code
# Likelihood Ratio Test
anova(exp_AFT, wei_AFT_surv)
  Terms Resid. Df    -2*LL Test Df Deviance Pr(>Chi)
1 treat        40 217.0481      NA       NA       NA
2 treat        39 213.1590    =  1 3.889116   0.0486

We could also use the lrtest() function from the lmtest() package to compare two nested flexsurvreg objects:

Show R code
# Make sure you have installed the lmtest package
# install.packages("lmtest")
# Load the lmtest package
library(lmtest)

# Likelihood Ratio Test

# It is recommended to list the most complex model first
lrtest(exp_PH, wei_PH)
Likelihood ratio test

Model 1: Surv(time, cens) ~ treat
Model 2: Surv(time, cens) ~ treat
  #Df  LogLik Df  Chisq Pr(>Chisq)  
1   2 -108.52                       
2   3 -106.58  1 3.8891     0.0486 *
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Also note that while the Wald Test (the p-value for Log(scale) in survreg output is a quick shortcut, as demonstrated above, the LRT is generally considered more robust for model selection.

2. Likelihood Ratio Test (LRT)

While the Likelihood Ratio Test (LRT) is useful for comparing nested models (like Exponential vs. Weibull), the Akaike Information Criterion (AIC) provides a more general framework for model selection. AIC allows us to compare any set of models, even those that are not nested, by penalizing for the number of parameters used.

While LRT asks: “Is the improvement in log-Likelihood large enough to be ‘real’?”, AIC asks: “Is the improvement in log-Likelihood large enough to justify the extra parameter we added?”

The AIC is defined as:

\[AIC = -2\log(L) + 2p,\]

where \(\log(L)\) is the log-likelihood of the model and \(p\) is the number of estimated parameters in the model.

The goal of AIC is to find the model that explains the most variation using the fewest possible parameters (the principle of parsimony). The \(+2p\) term penalises the model for adding extra parameters. This prevents overfitting. When comparing models, the one with the lower AIC value is considered the better fit.

While the flexsurvreg outputs for the Exponential and Weibull PH and AFT models automatically print out AIC, survreg output only prints out the log-likelihood values. However, we can extract the log-likelihood values from the survreg output and manually calculate the AICs to compare the two models:

Show R code
# Calculate manually

# Using Survreg output above that provides the log-likehood values
# The necessary values can be extracted in the following way:
# the first value of loglik element gives the log-likelihood for the model with only an intercept
# the second value gives the log-likelihood for the model with the covariates
logL_exp <- exp_AFT$loglik[2]
logL_exp
[1] -108.524
Show R code
logL_wei <- wei_AFT_surv$loglik[2]
logL_wei
[1] -106.5795
Show R code
# Calculate the AIC
AIC_exp <- -2 * (logL_exp) + 2*2
AIC_exp
[1] 221.0481
Show R code
AIC_wei <- -2 * (logL_wei) + 2*3
AIC_wei
[1] 219.159
Show R code
Delta_AIC <- AIC_exp - AIC_wei
Delta_AIC
[1] 1.889116

We can also use the AIC() function from the stats package:

Show R code
# Comparing AIC
AIC(exp_PH, wei_PH)
       df      AIC
exp_PH  2 221.0481
wei_PH  3 219.1590
Show R code
aic_table <- data.frame(
  Model = c("Exponential", "Weibull"),
  AIC = c(AIC(exp_PH), AIC(wei_PH))
)

# Sort by AIC (lowest is best)
aic_table <- aic_table[order(aic_table$AIC), ]
print(aic_table)
        Model      AIC
2     Weibull 219.1590
1 Exponential 221.0481

The Weibull model has a lower AIC (\(219.16 < 221.05\)), indicating that the improvement in log-likelihood gained by adding the shape parameter (\(\gamma\)) is large enough to justify the extra complexity. Thus, the Weibull model is the preferred fit for these data.

A ‘rule of thumb’ in biostatistics is that an AIC difference (\(\Delta AIC\)) of more than 2 is generally considered a meaningful improvement. In this data example, \(\Delta AIC = 1.89\), which is a borderline but justifiable reason to choose the Weibull model.

14 Log-logistic distribution

While the Weibull distribution is versatile, it is limited by the fact that its hazard function is monotonic — it must always increase, always decrease, or remain constant.

In many clinical scenarios, the risk of an event increases initially (e.g., following a major surgery such as transplantation) and then decreases as the patient recovers. For these situations, the Log-logistic distribution offers the necessary flexibility.

  • Hazard function:

\[h(t) = \frac {\lambda \gamma t^{\gamma-1}} {1+\lambda t^\gamma}\]

  • Survival function:

\[S(t) = \frac {1} {1+\lambda t^\gamma}\]

  • Probability density function:

\[f(t) = \frac {\lambda \gamma t^{\gamma-1}} {(1+\lambda t^\gamma)^2}\]

If \(\gamma \le 1\) hazard decreases over time and if \(\gamma > 1\) hazard first increases and then decreases over time:

Show R code
# Log-logistic hazard, survival and probability density functions
# assumes scale parameter lambda is 1

data_loglogistic <- data.frame(expand.grid(Time = seq(0.01, 3, 0.01),
                                           Gamma = c(0.5, 1, 2, 5))) %>%
  mutate(
    # Log-logistic Hazard
    Hazard = (Gamma*Time^(Gamma - 1)) / (1 + Time^Gamma),

    # Log-logistic Survival
    Survival = 1 / (1 + Time^Gamma),
    
    # Log-logistic PDF
    PDF = (Gamma*Time^(Gamma - 1)) / (1 + Time^Gamma)^2,
    
    c = recode_factor(Gamma, `0.5` = "gamma == 0.5", 
                      `1`   = "gamma == 1.0", 
                      `2`   = "gamma == 2.0", 
                      `5`   = "gamma == 5.0")
  )

# Plotting the hazard (the most interesting part of log-logistic)
ggplot(data_loglogistic, aes(x = Time, y = Hazard, color = c)) +
  geom_line(aes(linetype = c), linewidth = 0.8) +
  labs(title = "Log-logistic Hazard Functions",
       subtitle = "Notice the unimodal (up-then-down) shape when gamma > 1",
       color = "", linetype = "") +
  scale_colour_discrete(labels = scales::parse_format()) +
  scale_linetype_discrete(labels = scales::parse_format()) +
  theme_bw()

Show R code
# Plotting the survival
ggplot(data_loglogistic, aes(x = Time, y = Survival, color = c)) +
  geom_line(aes(linetype = c), linewidth = 0.8) +
  labs(title = "Log-logistic Survival Functions",
       color = "", linetype = "") +
  scale_colour_discrete(labels = scales::parse_format()) +
  scale_linetype_discrete(labels = scales::parse_format()) +
  theme_bw()

Show R code
# Plotting probability density functions
ggplot(data_loglogistic, aes(x = Time, y = PDF, color = c)) +
  geom_line(aes(linetype = c), linewidth = 0.8) +
  labs(title = "Log-logistic Probability Density Functions",
       color = "", linetype = "") +
  scale_colour_discrete(labels = scales::parse_format()) +
  scale_linetype_discrete(labels = scales::parse_format()) +
  theme_bw()

15 Data example: Fitting a log-logistic model

1. Log-logistic accelerated failure time (AFT) model

The log-logistic distribution accommodates an AFT model but not a PH model.

This can be seen by deriving the hazard ratio for this model:

\[ HR = \frac {\lambda_1} {1+\lambda_1 t^\gamma} \frac {1+\lambda_0 t^\gamma} {\lambda_0}, \] as the time variable \(t\) does not cancel out.

As the hazard function in log-logistic model has a complex (often unimodal) shape, the ratio of hazards between two groups typically changes as time progresses. The log-logistic model can be used when we suspect non-proportional hazards (e.g., when the treatment effect is very strong early on but fades later).

Fitting log-logistic AFT model

We can fit the log-logistic AFT model using flexsurvreg function from the flexsurv package, specifying dist="llogis":

Show R code
# LOG-LOGISTIC AFT MODEL

# Fitting log-logistic AFT model using flexsurvreg
llog_AFT <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "llogis")
# gives estimate for shape parameter and scale parameter
# gives alpha, TR and its 95% confidence interval
llog_AFT
Call:
flexsurvreg(formula = Surv(time, cens) ~ treat, data = leukdata, 
    dist = "llogis")

Estimates: 
           data mean  est    L95%   U95%   se     exp(est)  L95%   U95% 
shape         NA      1.830  1.363  2.456  0.275     NA        NA     NA
scale         NA      6.637  4.418  9.970  1.378     NA        NA     NA
treat6-MP  0.500      1.265  0.627  1.904  0.326  3.545     1.872  6.711

N = 42,  Events: 30,  Censored: 12
Total time at risk: 541
Log-likelihood = -107.6614, df = 3
AIC = 221.3228
Show R code
# This will calculate the actual hazard rates for each group at different times
summary(llog_AFT, newdata = data.frame(treat = c("control", "6-MP")), type = "hazard")
treat=control 
   time        est        lcl        ucl
1     1 0.05559348 0.01710429 0.12181596
2     2 0.09168701 0.03945990 0.16125169
3     3 0.11560823 0.05652333 0.19145226
4     4 0.12973570 0.07097029 0.21131321
5     5 0.13658720 0.08074732 0.21878217
6     6 0.13843226 0.08368432 0.22160157
7     7 0.13704689 0.08540003 0.21798510
8     8 0.13370164 0.08461501 0.21380256
9     9 0.12925486 0.08340716 0.20131671
10   10 0.12426427 0.08184300 0.18926043
11   11 0.11908127 0.07835245 0.17945021
12   12 0.11392066 0.07493404 0.17072982
13   13 0.10890879 0.07181310 0.16293036
14   15 0.09957569 0.06610480 0.14657526
15   16 0.09530154 0.06384999 0.13978512
16   17 0.09129243 0.06154118 0.13272110
17   19 0.08403094 0.05727156 0.12074770
18   20 0.08075081 0.05517645 0.11544348
19   22 0.07481425 0.05139265 0.10674939
20   23 0.07212723 0.04965875 0.10260388
21   25 0.06724484 0.04651627 0.09512522
22   32 0.05413207 0.03808962 0.07525846
23   34 0.05123425 0.03610421 0.07097718
24   35 0.04989393 0.03518540 0.06907393

treat=6-MP 
   time         est         lcl        ucl
1     1 0.005643366 0.001014487 0.02009116
2     2 0.009951271 0.002602659 0.02718494
3     3 0.013765985 0.004292289 0.03366425
4     4 0.017207751 0.005932107 0.03999970
5     5 0.020321906 0.007601904 0.04442905
6     6 0.023131853 0.009052069 0.04935033
7     7 0.025653632 0.010461474 0.05336424
8     8 0.027901194 0.011889943 0.05696108
9     9 0.029888446 0.013217420 0.06035663
10   10 0.031629968 0.014500510 0.06441320
11   11 0.033141120 0.015579707 0.06630824
12   12 0.034437872 0.016624292 0.06800508
13   13 0.035536527 0.017078561 0.06992289
14   15 0.037204614 0.018225786 0.07089317
15   16 0.037805703 0.018697303 0.07143751
16   17 0.038271586 0.018993311 0.07105677
17   19 0.038853148 0.020440149 0.07031863
18   20 0.038994214 0.020667774 0.06956018
19   22 0.039033113 0.021180538 0.06952122
20   23 0.038950545 0.021432724 0.06877540
21   25 0.038623595 0.022041138 0.06699623
22   32 0.036426708 0.021455874 0.06082470
23   34 0.035641995 0.021211027 0.05862058
24   35 0.035238294 0.021034054 0.05791538

We can also fit the log-logistic AFT model using survreg function from the survival package, specifying dist="loglogistic":

Show R code
# LOG-LOGISTIC AFT MODEL

# fitting log-logistic AFT model using survreg
llog_AFT_surv <- survreg(Surv(time, cens) ~ treat, data = leukdata, dist = "loglogistic")
summary(llog_AFT_surv)

Call:
survreg(formula = Surv(time, cens) ~ treat, data = leukdata, 
    dist = "loglogistic")
             Value Std. Error     z       p
(Intercept)  1.893      0.208  9.12 < 2e-16
treat6-MP    1.265      0.326  3.89   1e-04
Log(scale)  -0.604      0.150 -4.02 5.7e-05

Scale= 0.547 

Log logistic distribution
Loglik(model)= -107.7   Loglik(intercept only)= -115.4
    Chisq= 15.38 on 1 degrees of freedom, p= 8.8e-05 
Number of Newton-Raphson Iterations: 4 
n= 42 
Show R code
# Function predict can compute median survival time 
# Median survival time for placebo and treatment group
median_times <- predict(llog_AFT_surv, newdata = data.frame(treat = c("control", "6-MP")), 
                        type = "quantile", p = 0.5)
median_times
        1         2 
 6.637203 23.527124 

The shape parameter (\(\hat{\gamma} = 1.83\)) is reported on the natural scale in both outputs. As \(\gamma > 1\), the hazard is unimodal - it rises to a peak and then declines.

The log-logistic model output from flexsurvreg for the scale parameter (\(\hat{\gamma} = 6.637\)) gives the median survival time for the baseline group (placebo), i.e., 50% of the people in placebo group are expected to have experienced relapse by 6.6 weeks.

Consistent with the log-linear structure of the AFT model, the covariate effect \(\hat{\alpha_1} = 1.265\) is reported on the log-time scale. The flexsurvreg output also automatically provides the time ratio (TR)
\(TR = \exp(\hat{\alpha_1}) = \exp(1.265) = 3.55\) and its 95% confidence interval (1.87, 7.71).

Log-logistic proportional odds (PO) model

As demonstrated above, the log-logistic model does not generally satisfy the proportional hazards assumption because the hazard ratio varies over time. However, it satisfies the proportional odds (PO) assumption and is thus a proportional odds model.

A PO model is a model in which the survival odds is assumed to remain constant over time. In survival analysis, we can define odds similarly to how we do in logistic regression, as odds of an individual surviving beyond some time \(t\):

  • Survival Odds: \(\frac{S(t)}{1−S(t)}\)

  • Failure Odds: \(\frac{1-S(t)}{S(t)}\)

For the log-logistic distribution, where \(S(t) = \frac {1} {1+\lambda t^\gamma}\) and \(1-S(t) = \frac {\lambda t^\gamma} {1+\lambda t^\gamma}\), the survival odds is defined as:

\[\text{Odds}(t) = \frac {S(t)} {1−S(t)} = \frac {1} {\lambda t^\gamma} \] When we compare two groups (like treatment group 1 with \(\lambda_1\) and placebo group 0 with \(\lambda_0\)), the ratio of their survival odds looks like this:

\[\text{SOR} = \frac{\lambda_0 t^\gamma} {\lambda_1 t^\gamma} = \frac{\lambda_0} {\lambda_1} = \exp(-\beta_1),\] which is a constant as time \(t\) cancels out. This is why the proportional Odds assumption is satisfied for the log-logistic model.

We can estimate the SOR for our data example from the above output using the relationships \(\exp(-\beta_1) = \frac {\alpha_1} {\sigma} = \alpha_1 \gamma\), where the latter equation can be used for the flexsurvreg output and the either equation for the survreg output. We can also obtain SOR by writing out the SOR formula and applying it to the extended model output on survival estimates:

Show R code
# calculating SOR

# fitting log-logistic AFT model using flexsurvreg
# llog_AFT <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "llogis")

# fitting log-logistic AFT model using survreg
# llog_AFT_surv <- survreg(Surv(time, cens) ~ treat, data = leukdata, dist = "loglogistic")

# Extract the estimated coefficients
# 'res' contains the estimates in the 'est' column
estimates1 <- llog_AFT$res
#estimates1

estimates2 <- coef(llog_AFT_surv)
#estimates2

# Get the natural shape (gamma) and the AFT alpha from flexsurvreg output
gamma <- estimates1["shape", "est"] # flexsurvreg reports shape on natural scale
#gamma
alpha_aft <- estimates1["treat6-MP", "est"]
#alpha_aft

# Calculate the SOR
SOR1 <- exp(gamma * alpha_aft)
SOR1
[1] 10.12849
Show R code
SOR2 <- exp(llog_AFT_surv$coef[2] / llog_AFT_surv$scale)
SOR2
treat6-MP 
 10.12849 
Show R code
# This will calculate the survival probabilities for each group at different times
surv <- summary(llog_AFT, newdata = data.frame(treat = c("control", "6-MP")), type = "survival")

# this will calculate SOR using survival probability at time 1
# any time point can be used as SOR is a constant throughout the follow-up period
SOR3 = ( surv[[2]]$est[1] / (1-surv[[2]]$est[1]) ) / ( surv[[1]]$est[1] / (1-surv[[1]]$est[1]) )
SOR3
[1] 10.12849

Thus, assuming that log-logistic model fits the data and the proportional odds assumption is met, it appears that at any given time point, the odds of still being in remission for the treatment group are 10 times higher than the odds for the control (placebo) group.

15.1 Model validation

Again, before interpreting log-logistic model output, we must ensure the log-logistic distribution fits the data and that the accelerated failure time (AFT) assumption. and the proportional odds (PO) assumption are met.

We can carry out the same checks as before for Exponential and Weibull models:

1. Visual comparison with non-parametric estimates

First, we plot the fitted parametric survival curve \(\hat{S}(t)\) directly over the Kaplan-Meier (KM) curve.

Previously, we created custom plots based on the survival functions for the exponential and Weibull models. However, we can also use the flexsurvreg survival probability predictions in the following way:

Show R code
# 1. Fit the Kaplan-Meier (observed) by treatment group
km <- survfit(Surv(time, cens) ~ treat, data=leukdata)
km_df <- tidy(km)
# survival probability is now named estimate
#names(km_df)

# 2. Get the 'predicted' data (Log-Logistic)
# We generate predictions over a smooth sequence of time
time_seq <- seq(0, max(leukdata$time), length.out = 100)
pred_llogis <- summary(llog_AFT, t = time_seq, tidy = TRUE)
#pred_llogis

# 3. Create the Overlay Plot
ggplot() +
  # Observed: Kaplan-Meier steps
  geom_step(data = km_df, aes(x = time, y = estimate, color = strata), 
            linewidth = 0.8, alpha = 0.5) +
  # Predicted: Log-Logistic smooth lines
  geom_line(data = pred_llogis, aes(x = time, y = est, color = treat), 
            linewidth = 1.2) +
  labs(title = "Log-Logistic Model Fit: Predicted vs. Observed",
       subtitle = "Smooth lines = Log-Logistic Model; Stepped lines = Kaplan-Meier",
       x = "Time", y = "Survival Probability") +
  theme_minimal()

The log-logistic distribution appears to provide a reasonable good fit to the data.

We could also plot the hazard functions:

Show R code
plot(llog_AFT, type = "hazard", main = "Log-logistic hazard fit")

2. Linear transformation (graphical checks)

Second, we make comparisons with transformation of survival function, or in case the log-logistic model survival odds, into a linear form, and check whether the transformed data falls on a straight line.

For log-logistic model, the survival odds simplifies to

\[\text{Odds}(t) = \frac {S(t)} {1−S(t)} = \frac {1} {\lambda t^\gamma} = \lambda t^{-\gamma} \] and thus log(survival odds) \(= \log \lambda - \gamma \log t\). Therefore, the log survival odds, \(\log{\frac {S(t)} {1−S(t)}}\), is a linear function of log of time, with intercept \(\log{\lambda}\) and slope \(-\gamma\).

We can thus check whether plotting the log survival odds against the log of time gives a straight line of slope \(-\gamma\):

Show R code
library(tidyverse) # reloaded to be able to run this block independently
library(survival) # reloaded to be able to run this block independently

# Plotting survival odds against log of time

# 1. Fit the Kaplan-Meier (observed) by treatment group
km <- survfit( Surv(time, cens) ~ treat, data=leukdata)
km_df <- tidy(km)
# survival probability is now named estimate
#names(km_df)

# 2. Extract data and calculate log-odds
# Log-odds = log( S(t) / (1 - S(t)) )
km_data <- km_df %>%
  mutate(log_time = log(time),
      log_odds = log(estimate / (1-estimate))) %>%
  filter(is.finite(log_odds)) # Remove infinity values at S(t)=1 or 0

# 3. Create the plot
ggplot(km_data, aes(x = log_time, y = log_odds, color = factor(strata))) +
  geom_point() +
  geom_line() +
  geom_smooth(method = "lm", se = FALSE, linetype = "dashed") +
  labs(title = "Log-Logistic Diagnostic: Log-Odds Plot",
       subtitle = "Linearity = Log-logistic distribution; Parallelism = Proportional Odds",
       x = "log(Time)", 
       y = "log(S(t) / (1-S(t))",
       color = "Group") +
  theme_minimal()

Both plots look reasonably straight suggesting that the log-logistic assumption is reasonable.

We can also use this plot to assess proportional odds (PO) assumption. If the odds are proportional, the two lines should be parallel as they appear to be.

3. Graphical check for accelerated failure time assumption

We can evaluate whether the accelerated failure time (AFT) assumption for the log-logistic AFT model holds by checking whether there is a constant horizontal shift between survival curves plotted against log of time:

Show R code
library(broom) # reloaded to be able to run this block independently
library(tidyverse) # reloaded to be able to run this block independently
library(flexsurv) # reloaded to be able to run this block independently

# Plotting survival curves on log time scale to see the AFT model constant shift

# 1. Fit the Kaplan-Meier (observed) by treatment group
km <- survfit( Surv(time, cens) ~ treat, data=leukdata)
#summary(km)
#names(km)

# 2. Create a data frame for ggplot
plot_data <-broom::tidy(km)
#plot_data
# survival probability is now named estimate
#names(plot_data)

# 3. Fit log-logistic model
llog_AFT <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "llogis")
#exp_PH <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "exp")
#wei_AFT <- flexsurvreg(Surv(time,cens) ~ treat, data = leukdata, dist = "weibull")

# 4. Generate predicted survival probabilities for a range of times
# This creates a smooth curve for the log-time axis
time_seq <- seq(0.1, 35, length.out = 100)
pred_llogis <- summary(llog_AFT, t = time_seq, tidy = TRUE)
#pred_exp <- summary(exp_PH, t = time_seq, tidy = TRUE)
#pred_wei <- summary(wei_AFT, t = time_seq, tidy = TRUE)

# 5. Plot everything with a Log-transformed X-axis

# for log-logistic model only
ggplot() +
  # Kaplan-Meier steps (the raw data)
  geom_step(data = plot_data, aes(x = time, y = estimate, color = strata), linewidth = 0.8) +
  # Log-logistic fitted curves
  geom_line(data = pred_llogis, aes(x = time, y = est, color = treat), linetype = "dotdash", linewidth = 1) +
  scale_x_log10() + 
  labs(title = "KM vs. Log-logistic AFT fit on log time scale",
       x = "Log(Time)", y = "Survival Probability") +
  theme_bw()

Show R code
# also including exponential and Weibull models
ggplot() +
  # Kaplan-Meier steps (the raw data)
  geom_step(data = plot_data, aes(x = time, y = estimate, color = strata), linewidth = 0.8) +
  # Weibull fitted curves
  geom_line(data = pred_wei, aes(x = time, y = est, color = treat), linetype = "dotted", linewidth = 1) +
  # Exponential fitted curves
  geom_line(data = pred_exp, aes(x = time, y = est, color = treat), linetype = "dashed", linewidth = 1) +
  # Log-logistic fitted curves
  geom_line(data = pred_llogis, aes(x = time, y = est, color = treat), linetype = "dotdash", linewidth = 1) +
  scale_x_log10() + 
  labs(title = "KM vs. Parametric fits on log time scale",
       subtitle = "Dotted = Weibull, Dashed = Exponential, Dotdash = Log-logistic",
       x = "Log(Time)", y = "Survival Probability") +
  theme_bw()

The shift between survival curves does seem constant for log-logistic model.

15.2 Model selection

We can use AIC to compare the goodness-of-fit of exponential, Weibull and log-logistic models:

Show R code
aic_table <- data.frame(
  Model = c("Exponential", "Weibull", "Log-Logistic"),
  AIC = c(AIC(exp_PH), AIC(wei_AFT), AIC(llog_AFT))
)
# Sort by AIC (lowest is best)
aic_table <- aic_table[order(aic_table$AIC), ]
print(aic_table)
         Model      AIC
2      Weibull 219.1590
1  Exponential 221.0481
3 Log-Logistic 221.3228

Based on AIC, Weibull model and even exponential model appear to fit the data better than the log-logistic model although the differences in the AIC values between models are very small.

16 Log-normal distribution

The log-normal distribution is often used as an alternative to the log-logistic distribution because their hazard shapes are very similar. Both are flexible models within the accelerated failure time (AFT) framework.

A random variable \(T\) follows a log-normal distribution with parameters \(\mu\) (location) and \(\sigma\) (scale) if its natural logarithm, \(\log(T)\), follows a normal distribution with mean \(\mu\) and variance \(\sigma^2\):

\[\log(T) \sim N(\mu, \sigma^2)\] The probability density function of \(T\) is derived from the standard normal density \(\phi(z)\):

\[f(t) = \frac{1}{t\sigma} \phi\left( \frac{\log(t) - \mu}{\sigma} \right) = \frac{1}{t\sigma\sqrt{2\pi}} \exp\left[ -\frac{(\log(t) - \mu)^2}{2\sigma^2} \right]\] Unlike the exponential, Weibull and log-logistic models, the survival and hazard functions for the log-normal distribution cannot be expressed in simple algebraic forms. Instead, they are expressed in terms of the standard normal cumulative distribution function, \(\Phi(z)\):

\[S(t) = 1- \Phi\left( \frac{\log(t) - \mu}{\sigma} \right) = \int_t^\infty \frac{1}{x\sigma\sqrt{2\pi}} \exp\left[ -\frac{(\log(x) - \mu)^2}{2\sigma^2} \right]dx\]

\[h(t) = \frac{\frac{1}{t\sigma} \phi\left( \frac{\log(t) - \mu}{\sigma} \right)} {1- \Phi\left( \frac{\log(t) - \mu}{\sigma} \right)} = \frac{\frac{1}{t\sigma\sqrt{2\pi}} \exp\left[ -\frac{(\log(t) - \mu)^2}{2\sigma^2} \right]}{\int_t^\infty \frac{1}{x\sigma\sqrt{2\pi}} \exp\left[ -\frac{(\log(x) - \mu)^2}{2\sigma^2} \right]dx}\] Thus, hazard and survival functions can only be expressed in terms of integrals. Therefore, there is no closed-form solution for them. R software uses numerical approximation to calculate these values.

Like the log-logistic hazard where \(\gamma > 1\), the log-normal hazard is unimodal - it increases to a peak and then decreases towards zero:

Show R code
library(tidyverse) # reloaded to be able to run this code block independently

# Log-normal hazard, survival and probability density functions

# assumes mu = 0 (log1), meaning mean is t = 1

# Defining parameters: mu (meanlog) and sigma (sdlog)
# We will vary sigma to show how it affects the "peak" of the hazard
data_lognormal <- data.frame(expand.grid(Time = seq(0.01, 3, 0.01),
                                         Sigma = c(0.2, 0.5, 0.8, 1.2))) %>%
  mutate(
    # PDF: f(t)
    PDF = dlnorm(Time, meanlog = 0, sdlog = Sigma),
    
    # Survival: S(t) = 1 - CDF
    Survival = plnorm(Time, meanlog = 0, sdlog = Sigma, lower.tail = FALSE),
    
    # Hazard: h(t) = f(t) / S(t)
    Hazard = PDF / Survival,
    
    # Formatting labels for the plot legend
    s_label = recode_factor(Sigma, 
                            `0.2` = "sigma == 0.2", 
                            `0.5` = "sigma == 0.5", 
                            `0.8` = "sigma == 0.8", 
                            `1.2` = "sigma == 1.2")
  )

# 1. Plotting the Hazard Functions
ggplot(data_lognormal, aes(x = Time, y = Hazard, color = s_label)) +
  geom_line(aes(linetype = s_label), linewidth = 0.8) +
  labs(title = "Log-Normal Hazard Functions",
       subtitle = "The 'unimodal' shape: risk increases then eventually declines to zero",
       x = "Time (t)", y = "Hazard h(t)", color = "", linetype = "") +
  scale_colour_discrete(labels = scales::parse_format()) +
  scale_linetype_discrete(labels = scales::parse_format()) +
  theme_bw()

Show R code
# 2. Plotting the Survival Functions
ggplot(data_lognormal, aes(x = Time, y = Survival, color = s_label)) +
  geom_line(aes(linetype = s_label), linewidth = 0.8) +
  labs(title = "Log-Normal Survival Functions",
       x = "Time (t)", y = "Survival S(t)", color = "", linetype = "") +
  scale_colour_discrete(labels = scales::parse_format()) +
  scale_linetype_discrete(labels = scales::parse_format()) +
  theme_bw()

Show R code
# 3. Plotting the Probability Density Functions (PDF)
ggplot(data_lognormal, aes(x = Time, y = PDF, color = s_label)) +
  geom_line(aes(linetype = s_label), linewidth = 0.8) +
  labs(title = "Log-Normal Probability Density Functions",
       x = "Time (t)", y = "Density f(t)", color = "", linetype = "") +
  scale_colour_discrete(labels = scales::parse_format()) +
  scale_linetype_discrete(labels = scales::parse_format()) +
  theme_bw()

However, if the peak occurs very early (i.e., \(\sigma\) is large), the hazard may appear to be strictly decreasing over the actual range of the observed data:

Show R code
library(tidyverse) # reloaded to be able to run this code block independently

# Log-normal hazard for set of mu and sigma parameters, including large ones

params <- data.frame(
  mu =    c(0.1, 1, 2, 7),
  sigma = c(0.5, 0.5, 3, 5)
)

# Generate data with a Logarithmic sequence of Time to capture the peak
data_plot <- params %>%
  group_by(mu, sigma) %>%
  do(data.frame(Time = exp(seq(-5, 10, length.out = 1000)))) %>% 
  mutate(
    # PDF
    PDF = dlnorm(Time, meanlog = mu, sdlog = sigma),
    # Survival
    Survival = plnorm(Time, meanlog = mu, sdlog = sigma, lower.tail = FALSE),
    # Hazard: h(t) = f(t)/S(t)
    Hazard = PDF / Survival,
    Label = factor(paste0("mu=", mu, ", sigma=", sigma))
  )

# Plotting hazard functions
ggplot(data_plot, aes(x = Time, y = Hazard, color = Label)) +
  geom_line(linewidth = 1) +
  xlim(0, 10) + 
  labs(title = "Log-normal Hazard Functions",
       x = "Time (t)", y = "Hazard Rate h(t)") +
  theme_bw()

17 Data example: Fitting a log-normal model

The log-normal distribution accommodates an AFT model but not a PH model or a PO model. The covariates act as a multiplier on time, which matches the AFT assumption. However, When calculating the hazard ratio for a log-normal model, the \(t\) terms in the complex numerator and denominator do not cancel out. Consequently, the HR changes over time, violating PH assumption. Similarly, when calculating the survival odds ratio, the \(t\) terms also do not cancel out, violating the PO assumption.

Fitting log-normal AFT model

We can fit the log-normal AFT model using flexsurvreg function from the flexsurv package, specifying dist="lnorm":

Show R code
# LOG-NORMAL AFT MODEL

# Fitting Log-normal AFT model using flexsurvreg
ln_AFT <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "lnorm")
ln_AFT
Call:
flexsurvreg(formula = Surv(time, cens) ~ treat, data = leukdata, 
    dist = "lnorm")

Estimates: 
           data mean  est    L95%   U95%   se     exp(est)  L95%   U95% 
meanlog       NA      1.825  1.430  2.220  0.202     NA        NA     NA
sdlog         NA      0.924  0.713  1.197  0.122     NA        NA     NA
treat6-MP  0.500      1.347  0.727  1.967  0.316  3.845     2.068  7.150

N = 42,  Events: 30,  Censored: 12
Total time at risk: 541
Log-likelihood = -106.7046, df = 3
AIC = 219.4092
Show R code
# This will calculate the actual hazard rates for each group at different times
summary(ln_AFT, newdata = data.frame(treat = c("control", "6-MP")), type = "hazard")
treat=control 
   time        est        lcl       ucl
1     1 0.06286960 0.01296286 0.1505392
2     2 0.11455264 0.04806420 0.1891746
3     3 0.13473593 0.07317743 0.2059116
4     4 0.14127999 0.08538698 0.2135977
5     5 0.14189127 0.08839291 0.2193438
6     6 0.13981762 0.09000419 0.2147426
7     7 0.13652813 0.08859080 0.2115024
8     8 0.13272906 0.08571376 0.2058728
9     9 0.12877421 0.08296294 0.2012260
10   10 0.12484425 0.07989817 0.1955034
11   11 0.12103074 0.07666553 0.1896869
12   12 0.11737781 0.07443125 0.1855314
13   13 0.11390380 0.07217273 0.1813110
14   15 0.10750161 0.06719757 0.1732800
15   16 0.10456241 0.06556452 0.1690566
16   17 0.10178583 0.06390177 0.1651905
17   19 0.09667960 0.06086799 0.1575549
18   20 0.09432972 0.05917467 0.1541143
19   22 0.08998931 0.05616880 0.1476594
20   23 0.08798180 0.05513368 0.1446341
21   25 0.08425439 0.05278817 0.1389509
22   32 0.07361077 0.04610088 0.1223816
23   34 0.07110571 0.04464812 0.1188697
24   35 0.06992493 0.04396053 0.1169992

treat=6-MP 
   time         est          lcl        ucl
1     1 0.001190464 2.144008e-05 0.01128757
2     2 0.005924109 4.631242e-04 0.02368864
3     3 0.011747470 1.754281e-03 0.03357572
4     4 0.017129437 3.674298e-03 0.04031769
5     5 0.021645266 5.918740e-03 0.04609770
6     6 0.025287149 8.723714e-03 0.05070738
7     7 0.028170679 1.119757e-02 0.05440917
8     8 0.030431053 1.321153e-02 0.05702132
9     9 0.032190277 1.496967e-02 0.05800038
10   10 0.033549545 1.616983e-02 0.05931589
11   11 0.034590017 1.730337e-02 0.06071802
12   12 0.035375923 1.820430e-02 0.06143323
13   13 0.035957906 1.893328e-02 0.06216558
14   15 0.036661636 2.014988e-02 0.06285443
15   16 0.036840201 2.053926e-02 0.06314265
16   17 0.036931747 2.089471e-02 0.06428276
17   19 0.036915381 2.135138e-02 0.06408743
18   20 0.036831150 2.147665e-02 0.06360960
19   22 0.036554275 2.155799e-02 0.06283364
20   23 0.036374565 2.141731e-02 0.06261034
21   25 0.035956904 2.121064e-02 0.06244219
22   32 0.034187853 2.032770e-02 0.06015433
23   34 0.033651817 2.002525e-02 0.05929729
24   35 0.033383748 1.989617e-02 0.05873976

The parameters for a log-normal model in flexsurvreg are reported as follows:

  • meanlog (\(\hat{\mu}\)) is reported on the log-time scale. It represents the mean of the log-transformed survival times in the control (placebo) group.
  • sdlog (\(\hat{\sigma}\)) is reported on the natural scale to represent the actual standard deviation (\(\sigma\)) of survival times.
  • The coefficient (\(\hat{\alpha_1}\)) for the treatment covariate (treat6-MP) is
    reported on the log-time scale. As \(\hat{\alpha_1}>1\), treatment increases the survival time. To get the TR, we must again exponentiate the coefficient: \((\exp(1.347) = 3.85)\). The TR (exp(est) = 3.845) and its 95% confidence interval (2.068, 7.125) are automatically provided in the flexsurvreg output.

We can also fit the log-normal AFT model using survreg function from the survival package, specifying dist="lognormal":

Show R code
# LOG-NORMAL AFT MODEL

# Fitting Log-normal model using survreg
ln_AFT_surv <- survreg(Surv(time, cens) ~ treat, data = leukdata, dist = "lognormal")
summary(ln_AFT_surv)

Call:
survreg(formula = Surv(time, cens) ~ treat, data = leukdata, 
    dist = "lognormal")
              Value Std. Error     z       p
(Intercept)  1.8251     0.2016  9.05 < 2e-16
treat6-MP    1.3468     0.3165  4.26 2.1e-05
Log(scale)  -0.0792     0.1323 -0.60    0.55

Scale= 0.924 

Log Normal distribution
Loglik(model)= -106.7   Loglik(intercept only)= -115.4
    Chisq= 17.38 on 1 degrees of freedom, p= 3.1e-05 
Number of Newton-Raphson Iterations: 4 
n= 42 
Show R code
# Function predict can compute median survival time 
# Median survival time for placebo and treatment group
median_times <- predict(ln_AFT_surv, newdata = data.frame(treat = c("control", "6-MP")), 
                        type = "quantile", p = 0.5)
median_times
        1         2 
 6.203546 23.854377 

The survreg output is identical to the flexsurvreg output, with scale (\(\sigma\)) parameter reported both on the log-time scale and natural scale.

17.1 Model validation

Again, before interpreting log-normal model output, we must ensure the log-normal distribution fits the data and that the accelerated failure time (AFT) assumption. is met.

We can carry out similar checks as before for other models:

1. Visual comparison with non-parametric estimates

As before, we can plot the fitted parametric survival curve \(\hat{S}(t)\) directly over the Kaplan-Meier (KM) curve.

For a quick visual check, we can use plot(x) function. If x is a flexsurvreg object, R actually runs a function called plot.flexsurvreg(). It plots survival from a parametric model against non-parametric estimates, when you specify type = "survival":

Show R code
# Plot the Survival Function against the Kaplan-Meier data
plot(ln_AFT, type = "survival", main = "Log-Normal Survival Fit")

By plotting log-normal and log-logistic models in the same figure, we can see that they provide very similar fit:

Show R code
# Plot both to see which fits the data points better
plot(ln_AFT, type = "survival", ci = FALSE, col = "blue")
lines(llog_AFT, type = "survival", ci = FALSE, col = "red")
legend("topright", legend = c("Log-Normal", "Log-logistic"), col = c("blue", "red"), lty = 1)

Log-normal model appears to provide a reasonably good fit to the data.

We could also plot the hazard functions:

Show R code
plot(ln_AFT, type = "hazard", main = "Log-normal hazard fit")

2. Linear transformation (Graphical check)

If \(T\) follows a log-normal distribution, then \(\log(T)\) must follow a normal distribution. The Z score (normal quantile) is \(Z = \frac{\log{(t)}-\mu} {\sigma}\), which can be rearranged as \(Z = \frac{1}{\sigma} \log{(t)} - \frac{\mu} {\sigma}\).

We can thus check whether plotting the standard normal quantiles, calculated as \(\Phi^{-1}[1-S(t)]\), against \(\log(t)\), gives a straight line, with intercept \(-\frac{\mu} {\sigma}\) and slope \(\frac{1}{\sigma}\):

Show R code
library(tidyverse) # reloaded to be able to run this block independently
library(survival) # reloaded to be able to run this block independently

# Plotting normal quantiles against log of time 

# Kaplan-Meier fitted before: km

# 1. Fit the Kaplan-Meier (observed) by treatment group
km <- survfit(Surv(time, cens) ~ treat, data=leukdata)

# 2. Create a data frame for filtering and label the groups
# KM object often starts with a survival of 1 at time 0 and might end with a
# survival of 0
# gnorm(0) is -infinity and qnorm(1) is +infinity
# lm() below fails when it hits infinite values
# to fix this, we need to create a temporary data from, filter out those infinite
# points, and then plot the results
df <- data.frame(
  x=log(km$time), y=qnorm(1 - km$surv),
  # in survfit object, the information about the number of observations in each
  # group is stored in a vector called km$strata
  group = rep(names(km$strata), km$strata)
)

# 3. Remove the infinite values (where survival is 0 or 1)
df_clean <- df[is.finite(df$y), ]

# 4. Plot the cleaned data
plot(df_clean$x, df_clean$y, 
     main = "Log-Normal Q-Q Plot",
     xlab = "log(Time)", ylab = "Normal Quantiles"#,
     #pch = 19,
     #col ="blue"
     )

# 5. Extract and compare the slopes for each group
for(g in names(km$strata)){
  group_data <- subset(df_clean, group == g)
  fit_lm <- lm(y ~ x, data = group_data)
  
  coef(fit_lm)

# 6. Add a regression line for each group

  abline(fit_lm, col = ifelse(g == names(km$strata)[1], "red", "blue"))

# 7. Calculate sigma (1/slope) 
  slope <- coef(fit_lm)["x"]
  sigma_est <- 1 / slope
  print(paste("Group:", g, "| Estimated Sigma (1/slope):", round(sigma_est, 4)))
}
[1] "Group: treat=control | Estimated Sigma (1/slope): 1.0192"
[1] "Group: treat=6-MP | Estimated Sigma (1/slope): 1.3838"
Show R code
# 8. Add a legend
legend("topleft", legend = names(km$strata), col = c("red", "blue"), lty = 1, lwd = 2)

The plot of normal quantiles against \(\log(t)\) appear to produce relatively straight lines in both groups.

We can also use this plot to assess accelerated failure time assumption. If the AFT assumption is true, the treatment only changes the timing (\(\mu\)), not the spread or shape (\(\sigma\)). Therefore, the slopes for both groups must be equal.

Parallel lines on the plot are the visual proof that \(\sigma\) is constant across the groups. From the above plot it appears that the two lines are not exactly parallel. Above code also prints out the estimated \(\hat{\sigma}\) (1/slope) for both groups and they appear somewhat different, unlike the AFT model assumes. This implies that the treatment may be changing not only the timing but also the variance of the survival times. If this is the case, a single scale parameter \(\sigma\) is insufficient.

In flexsurvreg, it is possible to allow ancillary parameters like \(\sigma\) to depend on the treatment. You can do this using the anc argument. This tells R to estimate a separate \(\sigma\) for each group:

Show R code
# Fitting Log-normal model using flexsurvreg, 
# modelling both location and scale by treatment
ln_non_AFT <- flexsurvreg(Surv(time, cens) ~ treat, anc = list(sdlog = ~treat), data = leukdata, dist = "lnorm")
ln_non_AFT
Call:
flexsurvreg(formula = Surv(time, cens) ~ treat, anc = list(sdlog = ~treat), 
    data = leukdata, dist = "lnorm")

Estimates: 
                  data mean  est      L95%     U95%     se       exp(est)
meanlog                NA     1.8251   1.4388   2.2114   0.1971       NA 
sdlog                  NA     0.9032   0.6675   1.2222   0.1394       NA 
treat6-MP          0.5000     1.3779   0.6970   2.0589   0.3474   3.9667 
sdlog(treat6-MP)   0.5000     0.0803  -0.5056   0.6662   0.2989   1.0836 
                  L95%     U95%   
meanlog                NA       NA
sdlog                  NA       NA
treat6-MP          2.0076   7.8376
sdlog(treat6-MP)   0.6031   1.9468

N = 42,  Events: 30,  Censored: 12
Total time at risk: 541
Log-likelihood = -106.6677, df = 4
AIC = 221.3355
Show R code
# This will calculate the actual hazard rates for each group at different times
#summary(ln_non_AFT, newdata = data.frame(treat = c("control", "6-MP")), type = "hazard")

The output is otherwise similar to the AFT model output, but gives separate \(\sigma\) for control and treatment group, with the latter reported both on the log-time scale and natural scale (exp(est)). The two scale parameters for control group \(\hat{\sigma_0} = 0.9032\) and \(\hat{\sigma_1} = 1.0836\) do not appear very different and fit within each other’s 95% confidence intervals. However, we can carry out a formal LRT test to evaluate which model fits the data better.

3. Graphical check for accelerated failure time assumption

We can also evaluate whether the accelerated failure time (AFT) assumption for the log-normal AFT model holds by checking whether there is a constant horizontal shift between survival curves plotted against log of time:

Show R code
library(broom) # reloaded to be able to run this block independently
library(tidyverse) # reloaded to be able to run this block independently
library(flexsurv) # reloaded to be able to run this block independently

# Plotting survival curves on log time scale to see the AFT model constant shift

# 1. Fit the Kaplan-Meier (observed) by treatment group
km <- survfit( Surv(time, cens) ~ treat, data=leukdata)
#summary(km)
#names(km)

# 2. Create a data frame for ggplot
plot_data <-broom::tidy(km)
#plot_data
# survival probability is now named estimate
#names(plot_data)

# 3. Fit log-normal model (and other models for comparison)
#exp_PH <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "exp")
#wei_AFT <- flexsurvreg(Surv(time,cens) ~ treat, data = leukdata, dist = "weibull")
#llog_AFT <- flexsurvreg( Surv(time, cens) ~ treat, data = leukdata, dist = "llogis")
ln_AFT <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "lnorm")


# 4. Generate predicted survival probabilities for a range of times
# This creates a smooth curve for the log-time axis
time_seq <- seq(0.1, 35, length.out = 100)
#pred_exp <- summary(exp_PH, t = time_seq, tidy = TRUE)
#pred_wei <- summary(wei_AFT, t = time_seq, tidy = TRUE)
#pred_llogis <- summary(llog_AFT, t = time_seq, tidy = TRUE)
pred_lnorm <- summary(ln_AFT, t = time_seq, tidy = TRUE)

# 5. Plot everything with a Log-transformed X-axis

# for log-normal model only
ggplot() +
  # Kaplan-Meier steps (the raw data)
  geom_step(data = plot_data, aes(x = time, y = estimate, color = strata), linewidth = 0.8) +
  # Log-logistic fitted curves
  geom_line(data = pred_lnorm, aes(x = time, y = est, color = treat), linetype = "longdash", linewidth = 1) +
  scale_x_log10() + 
  labs(title = "KM vs. Log-normal AFT fit on log time scale",
       x = "Log(Time)", y = "Survival Probability") +
  theme_bw()

Show R code
# also including exponential, Weibull, and log-logistic models
ggplot() +
  # Kaplan-Meier steps (the raw data)
  geom_step(data = plot_data, aes(x = time, y = estimate, color = strata), linewidth = 0.8) +
   # Exponential fitted curves
  geom_line(data = pred_exp, aes(x = time, y = est, color = treat), linetype = "dashed", linewidth = 1) +
  # Weibull fitted curves
  geom_line(data = pred_wei, aes(x = time, y = est, color = treat), linetype = "dotted", linewidth = 1) +
  # Log-logistic fitted curves
  geom_line(data = pred_llogis, aes(x = time, y = est, color = treat), linetype = "dotdash", linewidth = 1) +
  # Log-normal fitted curves
  geom_line(data = pred_lnorm, aes(x = time, y = est, color = treat), linetype = "longdash", linewidth = 1) +
  scale_x_log10() + 
  labs(title = "KM vs. Parametric fits on log time scale",
       subtitle = "Dashed = Exponential, Dotted = Weibull, Dotdash = Log-logistic, Longdash = Log-normal",
       x = "Log(Time)", y = "Survival Probability") +
  theme_bw()

The shift between survival curves seems constant for log-normal model.

17.2 Model selection

Since the AFT model with common scale parameter (constrained: \(\sigma_0 = \sigma_1\)) for treatment and control group is nested within the non-AFT model with separate scale parameters (unconstrained: \(\sigma_0 \ne \sigma_1\)) for treatment and control groups, we can compare their fit using the likelihood ratio test:

Show R code
# 1. Fit the standard AFT model (common scale)
ln_AFT <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "lnorm")

# 2. Fit the extended non-AFT model (separate scales using 'anc')
ln_non_AFT <- flexsurvreg(Surv(time, cens) ~ treat, anc = list(sdlog = ~treat), data = leukdata, dist = "lnorm")

# 3. Perform the Likelihood Ratio Test
# Calculation: 2 * (LogLik(extended) - LogLik(standard))
lrt_stat <- 2 * (ln_non_AFT$loglik - ln_AFT$loglik)
p_val <- 1 - pchisq(lrt_stat, df = 1)

# 4. Display the results
cat("LRT Statistic:", round(lrt_stat, 4), "\n")
LRT Statistic: 0.0737 
Show R code
cat("P-value:", round(p_val, 4), "\n")
P-value: 0.786 

As the p-value is not significant, it appears that the slopes on the plot and in the non-AFT model were not different enough to require a more flexible model. We conclude that the treatment apparently influences only the timing (location) and not the variability (scale) of the event, supporting the AFT assumption, and not requiring us to model the ancillary parameters.

We can use AIC to compare the goodness-of-fit of exponential, Weibull, log-logistic, and log-normal models, including the AFT and non-AFT log-normal models:

Show R code
aic_table <- data.frame(
  Model = c("Exponential", "Weibull", "Log-Logistic", "Log-Normal AFT", "Log-Normal non-AFT"),
  AIC = c(AIC(exp_PH), AIC(wei_AFT), AIC(llog_AFT), AIC(ln_AFT), AIC(ln_non_AFT))
)
# Sort by AIC (lowest is best)
aic_table <- aic_table[order(aic_table$AIC), ]
print(aic_table)
               Model      AIC
2            Weibull 219.1590
4     Log-Normal AFT 219.4092
1        Exponential 221.0481
3       Log-Logistic 221.3228
5 Log-Normal non-AFT 221.3355

Based on AIC, the AFT log-normal model appears second best after the Weibull model and the non-AFT log-normal last, although the differences in the AIC values between models are again very small.

18 Gamma distribution

The shape of the gamma distribution is similar to that of the log-logistic or log-normal distributions, but they behave differently in their tails, as the tail of the gamma distribution decays more quickly.

If a random variable \(T\) follows a gamma distribution with parameters \(\lambda\) and \(k\), then

Probability density function:

\[f(t) = \frac{\lambda^k t^{k-1} e^{-\lambda t}}{\Gamma(k)},\] where the Gamma function is \(\Gamma(k) = \int_0^\infty x^{k-1} e^{-x} dx\).

Survival function:

\[S(t) = \frac{\Gamma(k, \lambda t)}{\Gamma(k)},\]

where \(\Gamma(k, \lambda t) = \int_{\lambda t}^\infty x^{k-1} e^{-x} dx\) is the upper incomplete Gamma function.

Hazard function:

\[h(t) = \frac{\lambda^k t^{k-1} e^{-\lambda t}}{\Gamma(k,\lambda t)},\] Hazard and survival functions cannot usually be expressed in simple closed forms and must be evaluated using integrals or numerical methods.

19 Generalised gamma distribution

The generalised gamma distributionThe generalised gamma distribution is a flexible three-parameter family (\(\lambda, k, p\)) that generalizes several key distributions.

If a random variable \(T\) follows a generalised gamma distribution with parameters \(\lambda\), \(k\), and \(p\) then

Probability density function:

\[f(t) = \frac{p \lambda^k t^{p k - 1} e^{-(\lambda t)^p}}{\Gamma(k)}\] Survival function:

\[S(t) = \frac{\Gamma(k, (\lambda t))^p}{\Gamma(k)},\]

where \(\Gamma(k, (\lambda t)^p = \int_{(\lambda t)^p}^\infty x^{k-1} e^{-x} dx\).

Hazard function:

\[h(t) = \frac{p \lambda^k t^{p k - 1} e^{-(\lambda t)^p}}{\Gamma(k, (\lambda t))^p}\] The power of this distribution lies in its nesting. It includes the following as special cases:

  • Exponential: when \(k=1\) and \(p=1\)
  • Weibull: when \(k=1\)
  • Lognormal: as \(k \to \infty\)
  • Standard Gamma: when \(p=1\).

20 Data example: Fitting gamma and generalised gamma models

Both the gamma and generalised gamma distributions are AFT models by nature.

Under specific conditions, they can accommodate PH models:

  • If \(k = 1\), the gamma model simplifies to the exponential model, which is a PH model.
  • If \(k = 1\) and \(p = 1\), the generalised gamma simplifies to the exponential model.
  • If \(k = 1\), the generalised gamma simplifies to the Weibull model, which is also a PH model.

These models are often used as non-proportional hazards models, allowing hazard ratios to change over time when the PH assumptions of simpler models fail.

We can plot the hazard functions, survival functions and probability density functions, showing exponential, Weibull and gamma models as special cases of the generalised gamma distribution:

Show R code
# Plotting hazard, survival and probability density functions

#library(flexsurv)
#library(ggplot2)
#library(tidyr)
#library(dplyr)

# 1. Define the parameter sets for comparison
# f(x; lambda, k, p)
# We will keep lambda (scale) constant at 1.0 for all
t <- seq(0.01, 5, length.out = 200)

# using flexsurv notation for generalised gamma model
models <- data.frame(
  label = c("Gen. Gamma (k=2, p=1.5)", "Gamma (p=1, k=2)", "Weibull (k=1, p=1.5)", "Exponential (k=1, p=1)"),
  # this is simplification that only works for Weibull and exponential
  mu = 0,        # mu = log(1/lambda). If lambda=1, mu=0. 
  sigma = c(1/(1.5*sqrt(2)), 1/sqrt(2), 1/1.5, 1), # sigma = 1/(p*sqrt(k))
  Q = c(1/sqrt(2), 1/sqrt(2), 1, 1)                # Q = 1/sqrt(k)
)

# 2. Generate the data for PDF, Survival, and Hazard
plot_data <- expand.grid(t = t, label = models$label) %>%
  left_join(models, by = "label") %>%
  mutate(
    PDF = dgengamma(t, mu = mu, sigma = sigma, Q = Q),
    Survival = pgengamma(t, mu = mu, sigma = sigma, Q = Q, lower.tail = FALSE),
    Hazard = hgengamma(t, mu = mu, sigma = sigma, Q = Q)
  )

# 3. Create the plots
# Function to generate uniform themed plots
make_plot <- function(data, y_var, title) {
  ggplot(data, aes(x = t, y = .data[[y_var]], color = label)) +
    geom_line(linewidth = 1.2) +
    labs(title = title, 
         x = "Time", 
         y = y_var, 
         color = "Distribution Family") +
    theme_minimal(base_size = 12) + # Slightly larger text for readability
    scale_color_brewer(palette = "Set1") #+
    #theme(legend.position = "bottom") # Moves legend to bottom for better horizontal space
}

p1 <- make_plot(plot_data, "Hazard", "Hazard Function (h(t))")
p2 <- make_plot(plot_data, "Survival", "Survival Function (S(t))")
p3 <- make_plot(plot_data, "PDF", "Probability Density Function (PDF)")

# Display plots
print(p1)

Show R code
print(p2)

Show R code
print(p3)

Fitting gamma AFT model

We can fit the gamma AFT model using flexsurvreg function from the flexsurv package, specifying dist="gamma":

Show R code
# GAMMA AFT MODEL

# Fitting gamma AFT model using flexsurvreg
g_AFT <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "gamma")
g_AFT
Call:
flexsurvreg(formula = Surv(time, cens) ~ treat, data = leukdata, 
    dist = "gamma")

Estimates: 
           data mean  est      L95%     U95%     se       exp(est)  L95%   
shape           NA     1.6539   1.0506   2.6035   0.3829       NA        NA
rate            NA     0.1908   0.1087   0.3350   0.0548       NA        NA
treat6-MP   0.5000    -1.2831  -1.8866  -0.6796   0.3079   0.2772    0.1516
           U95%   
shape           NA
rate            NA
treat6-MP   0.5068

N = 42,  Events: 30,  Censored: 12
Total time at risk: 541
Log-likelihood = -106.4417, df = 3
AIC = 218.8834
Show R code
# This will calculate the actual hazard rates for each group at different times
# summary(surv_g1, newdata = data.frame(treat = c("control", "6-MP")), type = "hazard")

# This gives you the Time Ratio
# untransformed coefficients
g_AFT$coefficients
     shape       rate  treat6-MP 
 0.5031215 -1.6562742 -1.2830901 
Show R code
exp(-g_AFT$coefficients[-c(1,2)])
treat6-MP 
 3.607771 

The parameters for a gamma model in flexsurvreg are as follows:

  • shape (\(\hat{k}\)): This determines the behavior of the hazard over time:
    • If shape > 1: The hazard is increasing
    • If shape < 1: The hazard is decreasing
    • If shape = 1: The hazard is constant (if the confidence interval for shape includes 1, model can likely be simplified to an exponential model)
  • rate (\(\hat{\lambda}\)): This is a scale parameter that relates to how stretched out the distribution is over time. The relationship between the rate \(\lambda\) and the location (\(\mu\)) in the AFT model is inverse. By default, flexsurv uses a log-link: (\(\mu\)) = \(-log(\lambda)\). Thus, a higher rate means that events happen faster (smaller \(\mu\)) and lower rate means that events happen slower (larger \(\mu\)). rate is thus the baseline log-rate for someone with all covariates zero.
  • The coefficient (\(\hat{\alpha_1}\)): Because flexsurv models the rate parameter (\(\lambda\)), the coefficients reported in the est column have a specific AFT meaning:
    • Positive coefficient: Increases the rate of failure (shorter survival time)
    • Negative coefficient: Decreases the rate of failure (longer survival time)

Therefore, in gamma model output exp(est) refers to acceleration factor (AF) rather than time ratio (TR). We can obtain \(TR = \frac{1}{\exp(\hat{\alpha_1})} = \exp(-(-1.2831)) = 3.60\).

Fitting generalised gamma AFT model

We can fit the generalised gamma AFT model using flexsurvreg function from the flexsurv package, specifying dist="gengamma":

Show R code
# GENERALISED GAMMA AFT MODEL

# Fitting generalised gamma AFT model using flexsurvreg
gg_AFT <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "gengamma")
gg_AFT
Call:
flexsurvreg(formula = Surv(time, cens) ~ treat, data = leukdata, 
    dist = "gengamma")

Estimates: 
           data mean  est     L95%    U95%    se      exp(est)  L95%    U95%  
mu             NA      2.063   1.386   2.741   0.346      NA        NA      NA
sigma          NA      0.830   0.540   1.278   0.183      NA        NA      NA
Q              NA      0.551  -0.807   1.910   0.693      NA        NA      NA
treat6-MP   0.500      1.306   0.680   1.932   0.319   3.691     1.973   6.903

N = 42,  Events: 30,  Censored: 12
Total time at risk: 541
Log-likelihood = -106.3895, df = 4
AIC = 220.779
Show R code
# This will calculate the actual hazard rates for each group at different times
# summary(gg_AFT, newdata = data.frame(treat = c("control", "6-MP")), type = "hazard")

This generalised gamma model output from flexsurvreg uses the so called Prentice parameterisation that relates it to the AFT model:

  • mu (\(\mu\)): This is the location parameter, that determines the ‘center’ of the survival time on a log scale for the reference (placebo group). The treatment covariate further shifts this (\(\mu\)).
  • sigma (\(\sigma\)): This is the scale parameter, that determines the spread or the variance of the log-survival times.
  • Q: This is the shape parameter, that controls the skewness and tail behavior. In a way, it thus handles the departure from log-normality. You can use it to evaluate whether a simpler two-parameter model might be justified:
    • Q = 1: Weibull (hazard changes monotonically, only up or only down)
    • Q = 0: Log-Normal (hazard rises to a peak and then declines)
    • Q > 0: Standard gamma (hazard typically increases and then levels off)
    • Q < 0: Inverse gamma

The parameterisation used before in our gamma formulas is called Stacy parameterisation and you can transform the above output into Stacy parameters using the following formulas:

Show R code
# Transforming flexsurv 'gengamma' output into stacy parameters used in the formulas

stacy_converter <- function(model) {
  # 1. Extract Prentice parameters from the model
  # We use the 'res' table which contains natural scale estimates
  mu    <- model$res["mu", "est"]
  sigma <- model$res["sigma", "est"]
  Q     <- model$res["Q", "est"]
  
  # 2. Apply the transformation math
  k_stacy <- 1 / (Q^2)
  p_stacy <- Q / sigma
  lambda_stacy <- exp(mu-(log(k_stacy)/p_stacy))
  
  # 3. Format into a nice table
  results <- data.frame(
    Parameter = c("k (Shape)", "p (Power)", "lambda (Rate/Scale)"),
    Stacy_Estimate = c(k_stacy, p_stacy, lambda_stacy),
    stringsAsFactors = FALSE
  )
  
  return(results)
}

# Usage:
 stacy_params <- stacy_converter(gg_AFT)
 print(stacy_params)
            Parameter Stacy_Estimate
1           k (Shape)      3.2889892
2           p (Power)      0.6640113
3 lambda (Rate/Scale)      1.3104385

However, you can also obtain the Stacy parameters by fitting the generalised gamma AFT model using flexsurvreg function from the flexsurv package, specifying dist="gengamma.orig":

Show R code
# GENERALISED GAMMA AFT MODEL

# Fitting gamma model using flexsurvreg
gg_AFT2 <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "gengamma.orig")
gg_AFT2
Call:
flexsurvreg(formula = Surv(time, cens) ~ treat, data = leukdata, 
    dist = "gengamma.orig")

Estimates: 
           data mean  est       L95%      U95%      se        exp(est)
shape            NA   6.66e-01  3.76e-02  1.18e+01  9.77e-01        NA
scale            NA   1.33e+00  2.36e-06  7.51e+05  9.01e+00        NA
k                NA   3.27e+00  2.13e-02  5.02e+02  8.39e+00        NA
treat6-MP  5.00e-01   1.31e+00  6.79e-01  1.93e+00  3.20e-01  3.69e+00
           L95%      U95%    
shape            NA        NA
scale            NA        NA
k                NA        NA
treat6-MP  1.97e+00  6.91e+00

N = 42,  Events: 30,  Censored: 12
Total time at risk: 541
Log-likelihood = -106.3895, df = 4
AIC = 220.779

In the Stacy version:

  • shape \((p)\): A power parameter that affects the spread and skewness
  • scale \((\lambda)\): Directly relates to the time scale
  • k \((k)\)): Affects the tail behavior and overall family of the distribution

Gamma or generalised gamma models cannot be fit using survreg.

20.1 Model validation

To ensure the gamma distribution fits the data and the accelerated failure time (AFT) assumption. is met, we carry out similar checks as before for other models:

1. Visual comparison with non-parametric estimates

We can plot the fitted parametric survival curve \(\hat{S}(t)\) directly over the Kaplan-Meier (KM) curve.

For a quick visual check, we can use plot(x) function. If x is a flexsurvreg object, R actually runs a function called plot.flexsurvreg(). It plots survival from a parametric model against non-parametric estimates, when you specify type = "survival":

Show R code
# Plot the gamma survival function against the Kaplan-Meier data
plot(g_AFT, type = "survival", main = "Standard Gamma Survival Fit")

Show R code
# Plot both the gamma and generalised gamma survival functions to see which fits 
# the data points better
plot(gg_AFT, type = "survival", ci = FALSE, col = "blue", main = "Model Comparison")
lines(g_AFT, type = "survival", ci = FALSE, col = "red")

# Add a legend
legend("topright", legend=c("Generalised Gamma", "Gamma"), 
       col=c("blue", "red"), lwd=2)

Show R code
# Plot additionally Weibul (best AIC so far) to see which fits the data points better
plot(gg_AFT, type = "survival", ci = FALSE, col = "blue", main = "Model Comparison")
lines(g_AFT, type = "survival", ci = FALSE, col = "red")
lines(wei_AFT, type = "survival", ci = FALSE, col = "green")

# 4. Add a legend
legend("topright", legend=c("Generalised Gamma", "Gamma", "Weibull"), 
       col=c("blue", "red", "green"), lwd=2)

We could also plot the hazard functions:

Show R code
plot(g_AFT, type = "hazard", main = "Gamma hazard fit")

Show R code
plot(gg_AFT, type = "hazard", main = "Generalised gamma hazard fit")

2. Linear transformation (Graphical check)

If \(T\) follows a gamma distribution with scale \(\lambda\) and shape \(k\), then \(\frac{T}{\lambda} \sim Gamma(k, \lambda = 1). The observed quantile of time\)t_p =Q(p;k,1)$, where \(t_p =\lambda * Q(p;k,1)\) is the theoretical quantile function for a standard gamma distribution

We can check whether plotting the observed survival times against theoretical gamma quantiles gives a straight line, with intercept 0 and slope :

Show R code
library(tidyverse) # reloaded to be able to run this block independently
library(survival) # reloaded to be able to run this block independently

# Plotting observed quantiles against theoretical gamma quantiles 

# 1. Fit your gamma model
g_AFT <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "gamma")

# 2. Extract the common shape parameter (k)
k_hat <- g_AFT$res["shape", "est"]
#rate_hat <- surv_g1$res["rate", "est"]
#a_hat <- 1 / rate_hat

# 3. Fit the Kaplan-Meier (observed) by treatment group
km <- survfit(Surv(time, cens) ~ treat, data=leukdata)

# 4. Create a data frame for filtering and label the groups
# KM object often starts with a survival of 1
# qgamma(1) is infinity
# lm() below fails when it hits infinite values
# to fix this, we need to create a temporary data from, filter out those infinite
# points, and then plot the results
df <- data.frame(
  time_obs = km$time, theo_quantiles = qgamma(1 - km$surv, shape = k_hat, scale = 1),
  # in survfit object, the information about the number of observations in each
  # group is stored in a vector called km$strata
  group = rep(names(km$strata), km$strata)
)

# 5. Remove the infinite values (where survival is 1)
df_clean <- df[is.finite(df$theo_quantiles), ]

# 6. Initialise the plot - created empty plot to add points in loop
plot(df_clean$theo_quantiles, df_clean$time_obs, type = "n",
xlim = c(0, max(df_clean$theo_quantiles)), 
ylim = c(0, max(df_clean$time_obs)),
xlab = "Theoretical Quantiles (Standard Gamma, Scale=1)",
ylab = "Observed survival time",
main = "Gamma Q-Q Plot",
xaxs = "i", yaxs = "i")

# 7. Extract slopes and add lines for each group
group_names <- names(km$strata)
cols <- c("blue", "red")

for(i in 1:length(group_names)){
  g <- group_names[i]
  group_data <- subset(df_clean, group == g)
  
   # Plot points for this group
  points(group_data$theo_quantiles, group_data$time_obs, col = cols[i], pch = 19)
  
  # Fit a regression line through origin (Time = scale * theoretical gamma quintile)
  # The slope of this line is the group-specific scale (a)
  # '-1' forces the intercept to be 0
  fit_lm <- lm(time_obs ~ theo_quantiles - 1, data = group_data)
  
  # Add the line
  abline(fit_lm, col = cols[i], lwd = 2)

# 7. Slope is the estimated scale a 
  scale_est <- coef(fit_lm)["theo_quantiles"]
  print(paste("Group:", g, "| Estimated scale:", round(scale_est, 4)))
}
[1] "Group: treat=control | Estimated scale: 4.9925"
[1] "Group: treat=6-MP | Estimated scale: 18.6808"
Show R code
# 8. Add a legend
legend("topleft", legend = group_names, col = cols, lty = 1, lwd = 2, pch = 19)

The plot of theoretical quantiles against observed survival times appear to produce relatively straight lines in both groups.

3. Graphical check for accelerated failure time assumption

We can also evaluate whether the accelerated failure time (AFT) assumption for the gamma and generalised gamma AFT models hold by checking whether there is a constant horizontal shift between survival curves plotted against log of time:

Show R code
library(broom) # reloaded to be able to run this block independently
library(tidyverse) # reloaded to be able to run this block independently
library(flexsurv) # reloaded to be able to run this block independently

# Plotting survival curves on log time scale to see the AFT model constant shift

# 1. Fit the Kaplan-Meier (observed) by treatment group
km <- survfit( Surv(time, cens) ~ treat, data=leukdata)
#summary(km)
#names(km)

# 2. Create a data frame for ggplot
plot_data <-broom::tidy(km)
#plot_data
# survival probability is now named estimate
#names(plot_data)

# 3. Fit log-normal model (and other models for comparison)
#exp_PH <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "exp")
#wei_AFT <- flexsurvreg(Surv(time,cens) ~ treat, data = leukdata, dist = "weibull")
#llog_AFT <- flexsurvreg( Surv(time, cens) ~ treat, data = leukdata, dist = "llogis")
#ln_AFT <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "lnorm")
g_AFT <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "gamma")
gg_AFT <- flexsurvreg(Surv(time, cens) ~ treat, data = leukdata, dist = "gengamma")


# 4. Generate predicted survival probabilities for a range of times
# This creates a smooth curve for the log-time axis
time_seq <- seq(0.1, 35, length.out = 100)
#pred_exp <- summary(exp_PH, t = time_seq, tidy = TRUE)
#pred_wei <- summary(wei_AFT, t = time_seq, tidy = TRUE)
#pred_llogis <- summary(llog_AFT, t = time_seq, tidy = TRUE)
#pred_lnorm <- summary(ln_AFT, t = time_seq, tidy = TRUE)
pred_g <- summary(g_AFT, t = time_seq, tidy = TRUE)
pred_gg <- summary(gg_AFT, t = time_seq, tidy = TRUE)

# 5. Plot everything with a Log-transformed X-axis

# for log-normal model only
ggplot() +
  # Kaplan-Meier steps (the raw data)
  geom_step(data = plot_data, aes(x = time, y = estimate, color = strata), linewidth = 0.8) +
  # Gamma fitted curves
  geom_line(data = pred_g, aes(x = time, y = est, color = treat), linetype = "dashed", linewidth = 1) +
  scale_x_log10() + 
  labs(title = "KM vs. Gamma AFT fit on log time scale",
       x = "Log(Time)", y = "Survival Probability") +
  theme_bw()

Show R code
# also including exponential, Weibull, and log-logistic models
ggplot() +
  # Kaplan-Meier steps (the raw data)
  geom_step(data = plot_data, aes(x = time, y = estimate, color = strata), linewidth = 0.8) +
   # Exponential fitted curves
  # geom_line(data = pred_exp, aes(x = time, y = est, color = treat), linetype = "dashed", linewidth = 1) +
  # Weibull fitted curves
  # geom_line(data = pred_wei, aes(x = time, y = est, color = treat), linetype = "dotted", linewidth = 1) +
  # Log-logistic fitted curves
  # geom_line(data = pred_llogis, aes(x = time, y = est, color = treat), linetype = "dotdash", linewidth = 1) +
  # Log-normal fitted curves
  # geom_line(data = pred_lnorm, aes(x = time, y = est, color = treat), linetype = "longdash", linewidth = 1) +
  # Gamma fitted curves
  geom_line(data = pred_g, aes(x = time, y = est, color = treat), linetype = "dashed", linewidth = 1) +
  # Generalised gamma fitted curves
  geom_line(data = pred_gg, aes(x = time, y = est, color = treat), linetype = "dotted", linewidth = 1) +
  scale_x_log10() + 
  labs(title = "KM vs. Parametric fits on log time scale",
       subtitle = "Dashed = Gamma, Dotted = Generalised Gamma", #Dotdash = Log-logistic, Longdash = Log-normal",
       x = "Log(Time)", y = "Survival Probability") +
  theme_bw()

The shift between survival curves seems constant for both models.

20.2 Model selection

In research you often don’t know which distribution fits your data best. Because generalised gamma distribution contains all these other distributions as special cases, you can:

  1. Fit a generalised gamma model
  2. Check the confidence intervals for \(k\) and \(p\)
  3. If p ≈ 1: You can simplify your model to a gamma model
  4. If k ≈ 1: You can simplify your model to a Weibull model
  5. If k ≈ 1, p ≈ 1: You can simplify your model to an exponential model

Because these other models are nested of with the generalised gamma model, we can carry out a formal Likelihood Ratio Test (LRT) between two models to determine which model provides a better fit for the data.

To test whether a generalised gamma \((\lambda, k, p\) distribution gives a significantly better fit than a Weibull \((\lambda, p\)) distribution, we test the following null hypothesis:

\[H_0: k = 1 \quad \text{(Weibull model)}\] \[H_1: k \neq 1 \quad \text{(Generalised gamma model)}\] Under the null hypothesis \(H_0\), the likelihood ratio test statistic \(W\) follows a chi-squared distribution with 1 degree of freedom (\(\chi^2_1\)):

\[W(\hat{\lambda}, \hat{p}, 1) = 2 \left[ \log L(\hat{\lambda}, \hat{p}, \hat{k}) - \log L(\hat{\lambda_0}, \hat{p_0}, 1) \right] \sim \chi_1^2,\] where \(\hat{\lambda}\), \(\hat{p}\), \(\hat{k}\) are the maximum likelihood estimates (MLEs) for a generalised gamma model with unconstrained \(k\), and \(\hat{\lambda}_0\), \(\hat{p}_0\) are the MLEs for a generalised gamma model with the constraint \(k = 1\) (i.e., MLE for the Weibull model).

The Weibull and generalised gamma model outputs we ran before provide the necessary log-likelihood values for carrying out the LRT. We can thus carry out the LRT for example in the following way:

Show R code
# Calculate manually

# Using flexsurvreg outputs that provide the log-likehood values
# The necessary values can be extracted in the following way:
logL_wei <- wei_AFT_surv$loglik
logL_wei
[1] -116.4054 -106.5795
Show R code
logL_gg <- gg_AFT$loglik
logL_gg
[1] -106.3895
Show R code
# Calculate the test statistic
LRT_statistic <- 2 * (logL_gg - logL_wei)
LRT_statistic
[1] 20.0318402  0.3800095
Show R code
# Calculate the P-value
# We use the Chi-squared distribution with 1 degree of freedom
p_value <- pchisq(LRT_statistic, df = 1, lower.tail = FALSE)
p_value
[1] 7.616337e-06 5.375982e-01

The non-significant p-value (\(p = 0.54\)) indicates that allowing the parameter \(k\) to vary (generalised gamma model) does not significantly improve the model’s log-likelihood compared to fixing it at 1 (Weibull model). We would thus fail to reject the null hypothesis at 5% significance level, and conclude that Weibull distribution is better than generalised gamma distribution in this case.

We could also use the lrtest() function from the lmtest() package to compare two nested flexsurvreg objects:

Show R code
# Make sure you have installed the lmtest package
# install.packages("lmtest")
# Load the lmtest package
library(lmtest)

# Likelihood Ratio Test

# It is recommended to list the most complex model first
lrtest(gg_AFT, wei_AFT)
Likelihood ratio test

Model 1: Surv(time, cens) ~ treat
Model 2: Surv(time, cens) ~ treat
  #Df  LogLik Df Chisq Pr(>Chisq)
1   4 -106.39                    
2   3 -106.58 -1  0.38     0.5376
Show R code
lrtest(gg_AFT, g_AFT)
Likelihood ratio test

Model 1: Surv(time, cens) ~ treat
Model 2: Surv(time, cens) ~ treat
  #Df  LogLik Df  Chisq Pr(>Chisq)
1   4 -106.39                     
2   3 -106.44 -1 0.1044     0.7466

Also when comparing the fit of standard gamma model and generalised gamma model, we fail to reject the null hypothesis, and also conclude that standard gamma distribution is better than generalised gamma distribution in this case.

We can use AIC to compare the goodness-of-fit of Weibull and standard gamma models, as well as other models we have fit:

Show R code
aic_table <- data.frame(
  Model = c("Exponential", "Weibull", "Log-Logistic", "Log-Normal",
            "Gamma", "Generalised Gamma"),
  AIC = c(AIC(exp_PH), AIC(wei_AFT), AIC(llog_AFT), AIC(ln_AFT),
          AIC(g_AFT), AIC(gg_AFT))
)
# Sort by AIC (lowest is best)
aic_table <- aic_table[order(aic_table$AIC), ]
print(aic_table)
              Model      AIC
5             Gamma 218.8834
2           Weibull 219.1590
4        Log-Normal 219.4092
6 Generalised Gamma 220.7790
1       Exponential 221.0481
3      Log-Logistic 221.3228

Based on AIC, the standard gamma model actually comes on top before Weibull, although the differences in AICs are very small.

We could also use BIC:

\(BIC = -2\log(L) + k \cdot \ln(n)\)

which uses a stricter penalty (whenever log(n) > 2, which happens once n > 7) and is less likely to ‘overfit’ by adding unnecessary parameters. This would only affect the ordering of models with different number of parameters (i.e., exponential and generalised gamm), not the models with the same number of parameters:

Show R code
bic_table <- data.frame(
  Model = c("Exponential", "Weibull", "Log-Logistic", "Log-Normal",
            "Gamma", "Generalised Gamma"),
  BIC = c(BIC(exp_PH), BIC(wei_AFT), BIC(llog_AFT), BIC(ln_AFT),
          BIC(g_AFT), BIC(gg_AFT))
)
# Sort by BIC (lowest is best)
bic_table <- bic_table[order(bic_table$BIC), ]
print(bic_table)
              Model      BIC
5             Gamma 224.0964
2           Weibull 224.3720
1       Exponential 224.5234
4        Log-Normal 224.6223
3      Log-Logistic 226.5358
6 Generalised Gamma 227.7297